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PREFACE 

This Second Annual Status Report describes the results of work performed dur- 
ing the first two years of the NASA Hot Section Technology program, "3-D 
Inelastic Analysis Methods for Hot Section Components" (contract NAS3-23697). 
The goal of the program is to develop computer codes which permit more accur- 
ate and efficient structural analyses of gas turbine blades, vanes, and com- 
bustor liners. The program is being conducted under the direction of Dr. C. C. 
Chamis of the NASA-Lewis Research Center. Prime contractor activities at 
United Technologies Corporation are managed by Dr. E. S. Todd. Subcontractor 
efforts at the United Technologies Research Center, MARC Analysis Research 
Corporation, and the State University of New York at Buffalo are led by Dr. B. 
N. Cassenti , Dr. J. C. Nagtegaal, and Dr. P. K. Banerjee, respectively. 
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SECTION 1.0 
INTRODUCTION 


Aircraft powerplant fuel consumption and expenditures for repair/replacement 
of worn or damaged parts make up a significant portion of commercial avia- 
tion's direct operating costs. For modern gas turbines, both factors depend 
heavily on the degree to which elevated flowpath temperatures are sustained in 
the hot section modules of the engine. Higher temperatures reduce fuel con- 
sumption by raising the basic efficiency of the gas generator thermodynamic 
cycle. At the same time, these elevated temperatures work to degrade the dura- 
bility of structural components (combustor liners, turbine blades and vanes, 
airseals, etc.) that must function adjacent to or within the hot gaspath it- 
self, leading in turn to larger maintenance/material costs. Pursuit of the 
best compromise between performance and durability presents a challenge that 
will continue to tax the ingenuity of advanced gas turbine design analysts for 
years to come. 

Hot section durability problems appear in a variety of forms, ranging from 
oxidation/corrosion, erosion, and distortion (creep deformations) to occur- 
rence of fatigue cracking. Even modest changes in shape, from erosion or dis- 
tortion of airfoils for example, can lead to measurable performance deteriora- 
tion that must be accurately predicted during propulsion system design to in- 
sure that long-term efficiency guarantees can be met. Larger distortions in- 
troduce serious problems such as hot spots and profile shifts resulting from 
diversion of cooling air, high vibratory stresses associated with loose tur- 
bine blade shrouds, difficult disassembly/reassembly of mating parts at over- 
haul, etc. These problems must be considered and efforts made to eliminate 
their effects during the engine design/development process. Initiation and 
propagation of fatigue cracks represents a direct threat to component struc- 
tural integrity and must be thoroughly understood and accurately predicted to 
insure continued safe and efficient engine operation. 
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Accurate prediction of component fatigue lives is strongly dependent on the 
success with which inelastic stress/strain states in the vicinity of holes, 
fillets, welds, and other discontinuities can be calculated. Stress/strain 
computations for hot section components are made particularly difficult by two 
factors - the high degree of geometrical irregularity which accompanies so- 
phisticated cooling schemes, and complex nonlinear material behavior associ- 
ated with high temperature creep/plasticity effects. Since cooling air ex- 
traction reduces engine cycle efficiency, concerted efforts are made to mini- 
mize its use with the result that elaborate internal passages and surface 
ports are employed to selectively bathe local regions (airfoil leading edges, 
louver liner lips, etc.) for which the high temperature environment is most 
severe. These cooling features frequently interrupt load paths and introduce 
complex temperature gradients to the extent that the basic assumptions of one- 
and two-dimensional stress analysis procedures are seriously compromised and 
the use of three-dimensional techniques becomes mandatory. Even in the pres- 
ence of cooling, component temperature and stress levels remain high relative 
to the material's melting point and yield strength values. The combinations of 
centrifugal, aerodynamic, thermal, and other mechanical loadings that typical- 
ly occur in flight operation then serve to drive the underlying material re- 
sponse beyond accepted limits for linear elastic behavior and into the regime 
characterized by inelastic, time-dependent structural deformations. Thus, an 
ability to account for both complexities, three-dimensional and inelastic 
effects, becomes essential to the design of durable hot section components. 


General purpose finite element computer codes containing a variety of three- 
dimensional (brick) elements and inelastic material models have been available 
for more than a decade. Incorporation of such codes into the hot section de- 
sign process has been severely limited by high costs associated with the ex- 
tensive labor/computer/time resources required to obtain reasonably detailed 
results. Geometric modeling systems and automated input/output data processing 
packages have received first attention from software developers in recent 
years and will soon mature to the point that previous over-riding manpower 
concerns will be alleviated. Prohibitive amounts of Central Processing Unit 
(CPU) time are still required for execution of even modest-size three-dimen- 
sional inelastic stress analyses, however, and is chief among the obstacles 
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remaining to be remedied. With today's computers and solution algorithms, mod- 
sis described by a few hundred displacement degrees of freedom commonly con- 
sume one to three hours of mainframe CPU time during simulation of a single 
thermomechanical loading cycle. A sequence of many such cycles may, of course, 
be needed to reach the stabilized conditions of interest. Since accurate 
idealizations of components with only a few geometrical discontinuities can 
easily contain several thousand degrees of freedom, inelastic analysis of hot 
section hardware with existing codes falls outside the realm of practicality. 

The Inelastic Methods Program addresses the need to develop more efficient and 
accurate three-dimensional inelastic structural analysis procedures for gas 
turbine hot section components. A series of new, increasingly rigorous, stand- 
alone computer codes is being created for the comprehensive numerical analysis 
of combustor liners, turbine blades and vanes. Theoretical foundations for the 
codes feature mechanics of materials models, special finite element models, 
and boundary element models. Heavy attention will be given to evolution of 
novel modeling methods that permit non-burdensome yet accurate representations 
of geometrical discontinuities such as cooling holes and coating cracks. A 
selection of constitutive relations has been provided for economical or so- 
phisticated description of inelastic material behavior as desired. Finally, 
advantages which accrue from application of the improved codes to actual com- 
ponents will be demonstrated by execution of benchmark analyses for which 
experimental data exist. 
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SECTION 2.0 


SUMMARY 

The 3-D Inelastic Analysis Methods program is divided into two 24-month seg- 
ments: a base program, and an option program to be exercised at the discretion 
of the Government. During the base program, a series of new computer codes em- 
bodying a progression of mathematical models (mechanics of materials, special 
finite element, boundary element) is being developed for the streamlined 
analysis of combustor liners, turbine blades and turbine vanes. These models 
will address the effects of high temperatures and thermal /mechanical loadings 
on the local (stress/strain) and global (dynamics, buckling) structural behav- 
ior of the three selected components. 

The first year (Task I) of the base program dealt with "linear" theory in the 
sense that stresses/strains and temperatures in generic modeling regions are 
linear functions of the spatial coordinates, and solution increments for load, 
temperature and/or time are extrapolated linearly from previous information. 
Three linear formulation computer codes, hereafter referred to as MOMM (Me- 
chanics of Materials Model), MHOST (MARC -HOST ), and BEST (Boundary Element 
Stress Technology), have been created and are described in more detail in the 
First Annual Status Report (NASA CR-1 74700). 

The second half of the base program (Task II), as well as the option program 
(Tasks IV and V), will extend the models to include higher-order representa- 
tions of deformations and loads in space and time and deal more effectively 
with collections of discontinuities such as cooling holes and coating cracks. 
Work on Task II (polynomial theory) has been completed, and the results are 
given in the third section of this Second Annual Status Report. 

2.1 CONSTITUTIVE MODELS 

Three increasingly rigorous constitutive relationships are employed by MOMM, 
MHOST, and BEST to account for nonlinear material behavior (creep/plasticity 
effects) in the elevated temperature regime. The simplified model assumes a 
bilinear approximation of stress-strain response and generally glosses over 
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the complications associated with strain rate effects, etc. (Section 3.1.1). 
The state-of-the-art model partitions time- independent (plasticity) and time- 
dependent (creep) phenomena in the conventional way, invoking the Mises yield 
criterion and standard (isotropic, kinematic, combined) hardening rules for 
the former and a power law for the latter (Section 3.1.2). Walker's viscoplas- 
tic theory, which accounts for the interaction between creep and plasticity 
that occurs under cyclic loading conditions, has been adopted as the advanced 
constitutive model (Section 3.1.3). 

2.2 MECHANICS OF MATERIALS MODEL 

In essence, the Mechanics of Materials Model (MOMM) is a stiffness method fi- 
nite element code that utilizes one-, two- and three-dimensional arrays of 
beam elements to characterize hot section component behavior. Limitations of 
such beam model representations are recognized, of course, but are fully ac- 
ceptable in view of the benefits of having a fast, easy to use, computation- 
ally efficient tool available for application during the early phases of com- 
ponent design. The full complement of structural analysis types (static, buck- 
ling, vibration, dynamics) is provided by MOMM, in conjunction with the three 
constitutive models mentioned above. Capabilities of the code have been tested 
for a variety of relatively simple problem discretizations (examples are pro- 
vided in Section 3.2.2). 

2.3 SPECIAL FINITE ELEMENT MODEL 

The MHOST (MARC -HOST ) code employs both shell and solid (brick) elements in a 
mixed method framework to provide comprehensive capabilities for investigating 
local (stress/strain) and global (vibration, buckling) behavior of hot section 
components. Development of the code has taken full advantage of the wealth of 
technical expertise accumulated at the MARC Corporation over the last decade 
in support of their own commercially available software packages to create 
new/improved algorithms (Section 3.3.4) that promise to significantly reduce 
CPU (central processing unit) time requirements for three-dimensional analy- 
ses. Second generation (Task II) MHOST code is operational and has been tested 
with a variety of academic as well as engine-related configurations (Section 
3.3.6). 
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2.4 ADVANCED FORMULATION (BOUNDARY ELEMENT) MODEL 


Development of the new BEST (Boundary Element stress Technology) co de consti- 
tutes a very important accomplishment of the Task II effort. The difficult 
challenge of extending the basic theory and algorithms to deal effectively 
with inelastic and dynamic effects in three-space was successfully met by com- 
bining the special skills and efforts of the research and programming teams at 
SUNY-B and P&W. As with MOMM and MHOST, the second version of BEST is opera- 
tional and has been exercised with a number of small and large test cases 
(Section 3.4.5). While MHOST and BEST are currently viewed as mutually comple- 
mentary, they are also competitors; and overall performance on large inelastic 
models will be watched with high interest as the codes continue to mature. 
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SECTION 3.0 
TECHNICAL PROGRESS 


3.1 CONSTITUTIVE MODELS 

Three material models are available for use with the mechanics of materials, 
special finite element, and boundary element models: 1) a simplified material 
model, 2) a state-of-the-art material model, and 3) an advanced material mod- 
el. The simplified model uses secant moduli and assumes a bilinear stress- 
strain response which is currently neither strain-rate nor temperature depen- 
dent. Later versions of the simplified material model may include provisions 
for both temperature and strain-rate dependence. The state-of-the-art material 
model is a standard elastic-plastic-creep model (Reference 1). The advanced 
model is a modified form of Walker's viscoplastic material model (References 2 
and 3). The following sections provide a detailed discussion of each of these 
model s . 

3.1.1 Simplified Secant Elastic Model 

In the simplified elastic model, stress-strain curves for various strain rates 
are the basic input material properties. Tension response is assumed to be the 
same as compression response. The initial response is represented by an elas- 
tic material with modulus, E Q , and Poisson's ratio, v Q . At the conclusion 
of the calculation for the response, an equivalent strain is predicted. At 
this strain, two equivalent stresses can be considered: 1) the calculated 
stress, and 2) the stress from the input stress- strain curves at the predicted 
strain. If the two stresses are sufficiently close in value, then the calcula- 
tions can be terminated. If the two stresses are not sufficiently close, then 
the new modulus is taken to be the stress from the stress-strain curves divid- 
ed by the strain, and the calculations are repeated. 


7 


This concept must now be expanded to multidimensional stress states. For this 
purpose, consider an elastic material, then: 


1 + V „ V _ x 

E ij F kk 6 ij 


(3.1-1) 


where: 


e .. is the mechanical strain tensor (i.e., total strain minus thermal strain), 

a-j j is the stress, and 

is the Kronecker delta. 

The stress and strain can be partitioned into deviatoric and volumetric parts, 

6j . = ♦ 1/3 t|(k «,j (3.1-2 

= S 1j * 1/3 ° kk 13-1 3 

The volumetric components, from equation (3.1-1) are related by: 


1 - 2v 


J kk = 3F °kk 


where K is the bulk modulus. 

The deviatoric parts can be shown to be related by: 

1 + V o 

e ij = E! S ij 

Let the equivalent stress be represented by: 

5 = V5T 2 = 73/2 sTT S.~ 

where J 2 1S the second invariant of the deviatoric stress tensor. 
Then, from equations (3.1-5) and (3.1-6): 

o = 1 + v >/3/2 e^ e^j 

Similarly, the equivalent strain can be taken to be: 

where j 2 is the second invariant of the deviatoric strain tensor. 
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Equation (3.1-7) now becomes: 


o * 



= 2Ge 


(3.1-9) 


Since only the ratio o/F will be used to represent the material response, an 
additional assumption is needed to obtain the second elastic constant. For 
this purpose, assume the bulk modulus is constant, and given by equation 
(3.1-4): 


K 


T(1 - 2v) = TTT 




2 (1 * v 0 ) G 0 
3 (1 - 2v 0 ) 


(3.1-10) 


where G Q , E Q , v Q are the moduli and Poisson's ratio at the origin (i.e., 
F=F*0 ) . The current shear modulus is known from the slope F/F. Then from equa- 
tion (3.1-9) 


2 (1 + v) G 

3 (1 - zvr 


Solving equation (3.1-11) for: 


2 (1 + v ) G 
o o 

3 (1 - 2v q ) 


V 



3 (G/G ) 

l 'Tn*^T7 

1 - 2v G 

o o 


(3.1-11) 


(3.1-12) 


figure 3.1-1 presents the variation in Poisson's ratio with modulus. The 
Young's modulus can be determined from: 


E = 2 (1 + v) G * (1 + v) a/l 
As an example, consider a uniaxial stress state: 


l a i=j=l 
ij 1 1 ° iVl. 


U 
Then: 


and 


( e i=J= 
■ = 7 -ve i=j= 

J ( o i ^ j 


E 1=j=l 

2,3 


(3.1-13) 


(3.1-14) 


o = a and e = (1 + v) e = e^j - e £2 


(3.1-15) 
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G/GO 

figure 3.1-1 Variation in Nu 
The "a, r curve is now the input stress-strain curve. 

To illustrate the convergence of the iterative procedure, consider three par- 
allel bars supporting an equivalent total load. The bars are assumed to be 
elastic-plastic. Each has a Young's modulus of 10 x 10 psi and a hardening 
slope of 0.5 x 10 6 psi. The yield stresses are different. The central bar 
will be assumed to have a yield stress of 20 ksi while the two outer bars have 
a yield stress of 10 ksi. The area of each bar is 1/3 in , making a total 
area of 1.0 in 2 . Figure 3.1-2 illustrates that convergence has occurred in 
six iterations for a total load of 30,000 lb and that each of the bars has 
yiel ded. 

The material constants for the simplified model are input to the computer code 
through data input cards. 
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figure 3.1-2 Three Bar Convergence Using Simplified Model 
3.1.2 Current State- of- the-Art Model 

The current state-of-the-art model has been taken to be the classical elastic- 
plastic-creep model that is available in the MARC code, and described in Ref- 
erence 1. The creep model is essentially a steady state power law (stress) 
model. The plasticity model includes isotropic, kinematic, and a combined 
hardening law. Both the creep and plasticity models assume no permanent volu- 
metric deformations. For the mechanics of materials computer code, the materi- 
al properties for the state-of-the-art constitutive model are included in data 
statements in subroutine SOACON. 

Plastic Iteration Procedure 


Consider the case of a small strain elastic-plastic response of a typical 
structure. Sufficiently large applied loads will result in permanent or plas- 
tic deformation. A procedure for calculating the response of the structure un- 
dergoing plastic deformation is required. 
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To evaluate the response of the structure, the loading history is divided into 
a number of incrementally applied loading steps. Each of these load increments 
can then be applied sequentially to the structure. An iterative scheme is then 
required to calculate the response of the structure to each individual load 
increment. 

At the beginning of a new load increment it may be assumed that the strain 
will change in a manner analogous to the previous increment. As an initial 
estimate all of the strain change is then assumed to be elastic. The change in 
the stresses can then be calculated using Hooke's Law or: 


A(J ij = L ijkl Ae kl 


( 3 . 1 - 16 ) 


where: ao-. is the incremental stress vector, 

Ae kl is the incremental total strain vector, and 


L 


e 

ijkl 


is the matrix of elastic constants. 


If the resulting total stress is within the yield surface, the matrix of mate- 
rial constants, L ij . kl is simply given by: 


L ijkl 



( 3 . 1 - 17 ) 


If the resulting total stress is outside the yield surface, weighted material 
constants and stiffness matrices will have to be calculated. It should be 
noted at this point that if a load increment is exceedingly large and if there 
is a sudden change in the type of loading, care must be taken in order to 
iterate to the correct solution. 
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If the resulting total stress is outside the yield surface, the fraction of 
the stress increment that remains elastic must be determined. This corresponds 
el 

to ae^I in Figure 3.1-3. If the yield surface in stress space is considered 
to be given by: 


then the appropriate m in: 

fio]: 1 + mao..) = 0 (3.1-18) 

I J I J 

may be determined where is the stress tensor from the previous increment. 

The mean material matrix is calculated from: 


L ijkl = mL i jkl + (1_m) L ijkl 
where is the tensor relating o.j and e kl . 


(3.1-19) 


Once the tensor j has been determined, standard solutions can be applied 
to find the incremental changes in the displacements, strains and loads. For 
example, if the strains are given by: 


|ae| = [B] { au| (3.1-20) 

where |au| is the vector of incremental nodal displacements, and [ 6 ] is the 
matrix relating the vector of element strains {ae} to the nodal displacements, 
the stiffness matrix can be found from: 


[K] = / Ce] T [D] [b] d V (3.1-21) 

v 

where [0] is the matrix representation of the tensor L.^^. The strain-dis- 
placement matrix [b] depends on the formulation of the element. 
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The Incremental nodal displacements and strains can be evaluated by solving 
for au in: 


[K] { au } = { ap} + |aG| (3.1-22) 

and then applying equation (3.1-20). 

In the mechanics of materials computer code the stiffness matrix K is held 
constant, and changes in the stiffness matrix are included in aG. 



Figure 3.1.3 Elastic-Plastic Strain Decompositions for Bilinear Stress-Strain 
Law 
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The term aP in equation (3.1-22) is the applied incremental load. The term aG 
is defined as the pseudo-load correction to the stiffness matrix due to in- 
elastic strains which is added to equation (3.1-22). The aG vector calculated 
from creep strain, for example, is shown in equation (3.1-34). 

One iteration cycle is completed each time the stiffness matrix is formed and 
the resulting equations solved. At the end of each cycle the resulting solu- 
tion must be tested for convergence. This is accomplished, by considering the 
change in energy. 


r 


e n _ e n-i e n _ e n-i 

? ’ 1/2 (E W ♦ E N_i ) 


(3.1-23) 


N-l 

where E is the change in energy summed over all elements on the previous 
N 

cycle and E is the energy including the present cycle. 


An accurate solution will usually result if r is maintained less than 0.1 for 
elastic-plastic problems. 

If the solution has satisfied the convergence, the stresses and strains can be 
updated and a new load increment added. If the solution has not converged, 
then a new guess for the strains, based on the latest cycle, must be input and 
the calculation procedure repeated. When the solution has not converged after 
a given number of cycles, the program should exit from the load incrementing 
loop. 


Figure 3.1-4 is a flow chart illustrating the small strain elastic-plastic it- 
eration procedure. 


For isotropic materials the moduli in equation (3.1-19) are given by: 


L ijkl = 1 + v { 6 ik 6 jl + 1 - V £v 6 ij 6 kll 


(3.1-24) 
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and 


L ijkl " T^T | 6 ik 5 jl + T - 27 6 ij 6 kl 

where E is Young's modulus 

v is Poisson's ratio 

6.. is the Kronecker delta 
^ J 

a. = H f + o w 


[ 


3/2 (S.j - ^•j j ) < S kl 
Z”'l + v 


o k i) 


1 + 2 


H + 


1 + v 

E 


>] 


r f p 

D • • = U £44 


n 


ij 


ij 


e P . are the plastic strain rates 
1 J 

t p L ,-> .p — 

e = ^2/3 

G is the kinematic hardening slope 
H is the isotropic hardening slope 


ay is the initial yield stress, and 
j = ai j - 1/3 aick«i j is the deviatoric stress. 


(3.1-25) 


(3.1-26) 

(3.1-27) 


(3.1-28) 


The strain rate has been decomposed into elastic (including thermal), plastic 
and creep components, or: 


i = e e + e P + e C (3.1-29) 

The plastic yield surface was assumed to satisfy an equivalent Mises yield 
surface given by: 


1/2 ( S ^ - Q i j ) (S.j - 8^) = 1/3 a 0 (3.1-30) 

The method presented in Reference 4 is used to calculate the elastic-plastic 
modul i . 
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ASSUME 
At =A< el + Af el Pl 


ASSUME 
Ac = A <7®^ — EA< 










Time Effects Iteration Procedure 


The creep strain rate will depend in general on the stress, the accumulated 
creep strain, the temperature and time. To illustrate the incrementing proce- 
dure, assume that the creep strain rate is normal to the Mises yield surface 
in stress space, then the creep strain rate is given by: 



p3/2 S„ S k1 j" 


3/2 S,j 


y/3/2~S, 


s 

mn mn 


(3.1-31) 


For a specific time increment the incremental creep strains were approximated 
by: 


. cr .cr A+ 
Ae ij = e ij At ‘ 

(3.1-32) 

The incremental displacements are: 


[K] )au( = {aP} + )aG( 

(3.1-33) 

where 


aG = /U] T [E] |Ae c } dV 

(3.1-34) 


is the pseudo-creep load, {a £ c } is the vector of element creep strain, and 
[E] is the elasticity matrix. The strain increment can be calculated from 
equation (3.1-20) and the strains, creep strains, stresses and displacements 
can be updated. 

A convergence test on the stresses should be performed. If the algorithm has 
not converged, a shorter time step should be used and the calculations re- 
peated. If the criterion has been satisfied, then the time step can be in- 
creased. figure 3.1-5 is a flow chart illustrating the small strain creep it- 
eration procedure. 
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3.1.3 Advanced Viscoplastic Model 


The viscoplastic model described in References 2 and 3 has been selected as 
the advanced constitutive model. Reference 2 describes the basic theory; while 
Reference 3 describes modifications to the form of the basic theory, and modi- 
fications to the material parameters for Hastelloy X. The modifications pro- 
vide more accuracy at relatively low temperatures. 

For uniaxial loading the viscoplastic material model (Figure 3.1-6) reduces to: 



sgn( o - o ) . -n) H - k > <°i> 

(3.1-35) 

o w - ko 


• • • 

£2 = n 2 c - n 3 |c| £2 

(3.1-36) 

C = e - £■ 

(3.1-37) 


where C is the inelastic strain, 

£2 is the back stress, 
o is the stress, 
e is the strain, and 

k, o n, K, n 2 , n 3 and E are material constants. 

The absolute value and unit ramp functions are represented by: 

i*i -{H $ 

and 

<x> - 1 0 x<0 

<x> " \ x x>_0 


(3.1-38) 

(3.1-39) 
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The inelastic strain in equation (3.1-35) consists of two components: 1) a 
time dependent power law creep component, containing the material constants n 
and k, and 2) a time independent plastic component, containing the material 
constants and k. The parameter a# becomes equivalent to the yield stress 
as: 1) k, in equation (3.1-35), approaches unity, and 2) the back stress, fi, 
approaches zero. The back stress is a key variable in many viscoplastic ma- 
terial models. Its evolution is given by equation (3.1-36). Equation (3.1-37) 
represents the inelastic strain as the difference between the total strain and 
the elastic strain. 


^ o (-f »ij"ni{)p-n)<g»jj«jj> 


H jj : (°i ♦ Cj; 


dn, . O / , dn, \ 0 

^'i* c 'i 7© 0 " (n ii“ n i,“ n,c i, , ( c_ "7 7e ®) * n '« 


K s K , — K j e s/3 


• f • D 

Cij 8 Cjj ♦ Cjj 




V ,S M X< kk* 2M ^r <r ij“ 8 li (3X * 2 '* ,o9,/2 '* 

-Tftjjftjj) 2 . 

«•/* =ij«!j • 

n . 3 5 [^i£tL . iiiii- - (iuSsiAl .( jSiliii- - g . 111 , 

^ C pq C pq C pq C PQ ' c pq c pq' ' C W ,C«»' * C 0a C O a ''*49 


•pq^pq 


*.j * 'ij ” J ®5j •’kk * 

iJ 2 s f(ivr a ii)(-Kr n M) 

o 

Material constants X t# » ,11 ,n,fn ,n, ,Oj , n s , n 4 , n 5 , n 6 , rt f , k , , k z , K , <r ^ depend on temperature 


Figure 3.1-6 Modified Walker's Theory 
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Subroutine HYPELA In the mechanics of materials computer code Integrates 
Walker's viscoplastic equations and calls subroutine HYPCON to evaluate the 
material parameters. HYPCON contains the latest estimates for the parameters 
in the modified Walker's theory. Each load increment in the analysis is divid- 
ed into NSPLIT subincrements. The integration of the constitutive equations is 
performed by using forward differences with a step size determined by dividing 
the load increment by NSPLIT. Subroutine HYPELA performs the integration in 
two ways: 1) a fixed step size, or 2) a variable step size. In the fixed step 
size, forward difference NSPLIT is the same for all load increments and sub- 
increments. 


In the variable step size, forward difference NSPLIT is determined by the mag- 
nitude of the change in a strain measure for every subincrement. The change in 
the strain measure is defined as: 


where 


V3ir 7 
E = aR + ~TT 


(3.1-40) 


A R = >/2/3 aC^ aC^ (3.1-41) 

aJ~ = 3/2 aS.. aS.. and (3.1-42) 

£ I J I J 


the quantity a T Is calculated and is stored as variable ERR0R0. There are 
three possible ways to determine NSPLIT. The method depends on the size of 
ERR0R0. If 


ERR0R2 < ERR0R0 < ERR0R1 , (3.1-43) 

then NSPLIT remains the same for the next subincrement (ERR0R1 and ERR0R2 are 
user-specified in HYPELA). If 


ERRORO < ERR0R2, 


(3.1-44) 


22 


4.u~ m ii cm tt -? _ j:...* j 

uicii iion.il 15 uiviucu in 

the nearest Integer. If 


two for the next subincrement and rounded (up) to 


ERRORO > ERR0R1, (3.1-45) 

then NSPLIT is doubled and the step is recomputed. The value of NSPLIT at the 
end of the increment is stored in the state variable TEMP (16). The initial 
value of NSPLIT is user-specified in HYPELA. The maximum value of NSPLIT is 
specified by MXSPLT. If NSPLIT exceeds MSXPLT, the message: 

"UNABLE TO REDUCE ERROR IN LESS THAN MXSPLT SUBINCREMENTS" 

is written where the value of MXSPLT is inserted in the WRITE statement. After 
this, the integration is performed using a constant step size. 
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3.1.4 List of Symbols 


List of Symbols 
Referenced Within Section 3.1 


Symbol 

Description 



Page 

e ij 

Strai n 



8 

°ij 

Stress 



8 

«1j 

Kronecker delta 



8 

V 

Poisson's ratio 



8 

E 

Young's modulus 



8 

e U 

Deviatoric strain 



8 

s ij 

Deviatoric stress 



8 

K 

Bulk modulus 



8 

J 2 

Second invariant of 
stress tensor 

the 

deviatoric 

8 

J2 

Second invariant of 
strain tensor 

the 

deviatoric 

8 


G 

L ijk£ 

{au| 

Cb] 

CK] 


Equivalent stress 

Equivalent strain 
Shear modulus 

Matrix of material constants 
Incremental nodal displacements 
Incremental stress 
Strain-displacement matrix 
Stiffness matrix 


8 

8 

9 

12 

13 

13 

13 

13 
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List of Symbols 
Referenced Within Section 3.1 


Symbol 

Description 

Page 

|ap| 

Incremental applied load vector 

14 

| aG^ 

Incremental pseudo-load vector 

14 

e n 

Energy in Nth cycle 

15 

r 

Convergence parameter 

15 

CE] 

Elasticity matrix 

18 

G 

Kinematic hardening slope 

16 

H 

Isotropic hardening slope 

16 

S2 

Back stress 

20 

c 

Inelastic strain 

20 

o 

x» ii i n, \ 



m, ni, n 2 , 1 

n3» n4, ns, \ 

Material constants 

20 


"O* "/» 'M* 1 

K2» k, / 
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3.2 MECHANICS OF MATERIALS MODEL 

3.2.1 Computer Program: Formulation/Description 

The three-dimensional nonlinear mechanics of materials finite element computer 
program utilizes an intersecting network of beams to model a structural com- 
ponent. The program calculates the total strain as a linear function of posi- 
tion in cross section and along the length of the beam. Three material models 
are included in the code: the simplified material model. Walker's viscoplastic 
material model, and the state-of-the-art material model. Static and transient 
analyses can be performed with applied loads, thermal loads, and enforced dis- 
placements. Frequencies and mode shapes using either initial or tangent stiff- 
ness is calculated; and buckling analysis is included in the static problem 
using initial or tangent stiffness. The program flow is summarized in Figure 
3.2-1. 

Input parameters to the computer code consist of information defining the 
model itself and information describing the method of solution desired. The 
model is defined by beams which are connected at grid points. The element co- 
ordinate system of a given beam is defined by an orientation grid point or 
vector. The geometry of a beam is rectangular in cross section, with the di- 
mensions of the cross section along the element coordinate axes specified. The 
material properties are specified for each beam, including Young's modulus, 
Poisson's ratio, mass density, coefficient of expansion, and yield stress. The 
initial temperature of the beam network is input, and the time at initial con- 
ditions is set to zero. A hardening slope for use with the simplified material 
model is entered, with a zero slope indicating perfectly-plastic behavior. 
Boundary conditions are specified by indicating at each node, a constrained or 
nonconstrained condition for the six degrees of freedom. 

Input associated with the selection of the method of solution include the pa- 
rameters that indicate: 

1. the choice of constitutive model to be used, 

2. the choice of a static or transient analysis. 
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3. the choice of Initial or tangent stiffness In solving for the fre- 
quencies and mode shapes, and 

4. the choice of Including buckling analysis with either Initial or tan- 
gent stiffness. 



Figure 3.2-1 3-D Inelastic Mechanics of Materials Computer Program Flow Chart 
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The number of integration points in each direction in each beam is user-speci- 
fied; stresses and strains are calculated at each integration point, and the 
user specifies up to 100 points along each element coordinate axis direction 
in each beam. The convergence value, defining the allowable energy change be- 
tween two consecutive iterations in the static analysis or allowable range in 
internal energy for the adaptive time step calculation in the transient analy- 
sis, is entered by the user. The number and type of loading increments are 
also specified. 

The stiffness and mass matrices for each beam in the element coordinate system 
are computed and transformed to the global coordinate system. The stiffness 
and mass matrices are then assembled to form the global mass and stiffness ma- 
trices. The boundary conditions are applied to the stiffness matrix, and the 
matrix is then inverted. Any change in the stiffness due to nonlinear effects 
will be accounted for in the pseudo-load vector; therefore, the stiffness ma- 
trix is only inverted once. 

Depending on user- input, the program now is directed to the appropriate branch 
of the program: static or transient analysis. For static analysis, the loading 
increment is read from the data input, including forces and moments or en- 
forced displacements, specified at each degree of freedom of the structure. 
The temperature increment is also entered. An initial incremental displacement 
vector is set to zero and strain, stress and pseudo-load vectors are calcu- 
lated from the incremental displacement vector using the mechanics of materi- 
als model selected by the user. The pseudo-load vector accounts for the 
effects of nonlinearity and allows the use of the original stiffness matrix 
throughout the calculations. The equations governing the system are as follows: 

[K] |au[ = |aP} + j aG } (3.2-1) 

where [K] = elastic stiffness matrix, 

{au| = incremental displacement vector, 

{ aP f = incremental applied load vector, and 

{aG} = pseudo-load vector, due to inelastic strains. 
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|aG| ■ / [B]^ [E] {ab} dV 


(3.2-2) 


where [B] = strain- displacement matrix, 

[E] = elasticity matrix, and 
|a 4 = Inelastic strains. 

Equation (3.2-1) is solved for the incremental displacement vector, au, which 
is substituted for the initial incremental displacement vector and used in the 
second iteration, continuing until the change in energy in two consecutive 
iterations is less than the convergence value input by the user. When conver- 
gence occurs, the incremental loading, displacements, strains and stresses for 
that loading increment are printed; the total load, displacement, strain and 
stress vectors, as well as temperature, are then updated. Each loading incre- 
ment is read in and executed similarly, and the values of stress, strain and 
displacement for the total loading are calculated and printed upon conclusion 
of the last increment. 

The transient analysis is based on a simple Euler integration and includes a 
self-adaptive time step scheme. Damping Is not Included directly in the tran- 
sient analysis but is present in the viscoplastic material models. The loading 
for each increment is the total load at that given time, which is entered into 
the program by a user-supplied subroutine. The temperature Increment and time 
step are also entered. As in the static branch, the initial displacement vec- 
tor is set to zero and the strains, stresses and pseudo-load vector are calcu- 
lated using the designated mechanics of material model. An Euler integration 
is then used to calculate current displacements at the end of the present time 
step. The governing equations are as follows: 


{Al = {af! - [K] |aU 0 [ 

(3.2-3) 

|aV} = [M]* 1 IaI * DT 

(3.2-4) 

{aU^ * ({v} ♦ 1/2 |aV|) * DT 

(3.2-5) 
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where |a| 

= acceleration vector. 

|aF[ 

= applied and pseudo-loads. 

CK] 

= elastic stiffness matrix. 

Ki 

= displacement vector at beginning of time step. 

|iV( 

= change- in- velocity vector. 

[M] 

= mass matrix. 

DT 

= time step. 

Kf 

= displacement vector at end of time step, and 

M 

= velocity vector. 


A measure of the work done and the change in internal energy of the system 
during the time step is computed, and the time step is adapted accordingly. If 
the time step is accepted, the current displacements, strains and stresses are 
printed, and the current displacements are inserted for the initial displace- 
ments in the following time step. If the time step is unacceptable according 
to the adaptive scheme, the time step is changed, the load is recalculated, 
and the displacements are reset to the initial value at the beginning of that 
time step. The analysis continues until the user-designated number of incre- 
ments are completed. 

Following the static or transient analysis, the user has a choice between 
calculating the lowest frequency and mode shape or all frequencies and mode 
shapes. The method of solution for the calculation of the lowest frequency and 
mode shape is the inverse power method, which is represented by the following 
expression: 


([K ] _1 [M] - X[I]) |x 1+1 [ = {x.| 


where [K] = stiffness matrix, 
[M] = mass matrix, 

[I] = identity matrix, 
X = eigenvalue, and 
| x [ = eigenvector. 


( 3 . 2 - 6 ) 


The method of solution in the calculation of all frequencies and mode shapes 
for a given problem is the Jacobi method, which is based on simple similarity 
transformations. 
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The procedure for determining the coefficients of the Inverse stiffness matrix 
Is one that can represent the original stiffness of the structure or the cur- 
rent stiffness Including nonlinear effects. A small load Is placed at one of 
the nonconstrained degrees of freedom of the structure, and the displacements 
are computed using the specified constitutive model. The coefficients of the 
appropriate row of the inverse stiffness matrix are calculated by dividing the 
calculated displacements by the applied force. This procedure is continued for 
each nonconstrained degree of freedom until an inverse stiffness matrix, with 
dimensions equal to the number of nonconstrained degrees of freedom of the 
structure, is formed. If the frequency is to be calculated using the initial 
stiffness of the structure, all variables used in the static or transient 
analysis are set to the original values. If the tangent stiffness is re- 
quested, all variables retain the current values for use in the frequency cal- 
culation. Only the initial stiffness option is available for use in a tran- 
sient analysis since current stiffness cannot be readily calculated. 

Buckling analysis can be executed in a static problem. The buckling analysis 
is based on a two step process similar to that in the NASTRAN finite element 
code. In the first step, the beam loads are determined. In the second step, a 
first order large displacement correction, proportional to the loads, is in- 
cluded in the stiffness matrix. Buckling occurs when the determinant of the 
new matrix vanishes. In the determination of the stiffness matrix used in the 
buckling calculation, the stiffness coefficients are calculated in the same 
fashion as was described in the frequency calculation, with the user choosing 
the initial or tangent stiffness. The beam loads are calculated using the 
initial stiffness matrix and then adding the pseudo-load vector. The actual 
buckling calculations are accomplished using the inverse power method to find 
the critical buckling factor and the buckled shape. 

Applied loads are entered by specifying the number of loading increments, the 
number of nodes with concentrated loads, and the number of beams with distrib- 
uted loads. Concentrated loads are input by defining the values of the six 
components of the loading at that node. Distributed loads are calculated using 
a consistent load formulation. The user inputs the six components of the load 
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at both end A and end 8 of the beam, and the distributed load is defined as 
linear between the two ends. A consistent load formulation is then used to 
calculate the nodal force vector for the beam, which is inserted into the 
appropriate section of the global force vector. Enforced displacements can be 
entered in the same way as concentrated loads; however, the corresponding 
boundary conditions must be constrained. 

Axial and transverse holes can be added to any beam. All holes are rectangular 
and the centerline of a hole must coincide with the centerline of the beam 
cross section, i.e., no off-center holes. Axial holes are input so that the 
wall thickness of the beam is constant. A section containing a hole is modeled 
as a beam with a constant cross section, and the values of cross sectional 
area, moments of inertia, and torsional stiffness are calculated and used to 
form the beam stiffness matrix in the same way as for a solid rectangular 
cross section. Calculations at integration points for determining the psuedo- 
load vector, strains, and stresses are bypassed if that integration point lies 
in a hole. 

Multipoint constraint equations can be supplied by the user. Each equation 
consists of a dependent degree of freedom defined in terms of the independent 
degrees of freedom. The equation is in the form: 

K\ - tV W (3 - 2 - 7) 

where ^ J = set of dependent degrees of freedom 
|u n J = set of independent degrees of freedom 
[G m ] = multipoint constraint matrix. 

The load- stiffness relationship can now be partitioned as follows: 


K nn K nm 
K mn 



( 3 . 2 - 8 ) 
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Writing equations (3.2-7) and (3.2-8) together and eliminating u gives: 


[K + K G m + G m T K 7, * G m T K Gl {uj = | P f + [G m T ] {Pj (3.2-9) 

nn nmm mnm m mm m ' n* ' n’ m ' m> 


or 


W * i p „i 


(3.2-10) 


where K = K 

nn 


nn 


K G 
nm m 


+ G T K T 
m nm 


+ G T K G 
m mm m 


= P + G T P 
nmm 


Cracks can be added to the structure. The user specifies the location and 
crack depth, and the computer code uses the elastic line spring model to de- 
termine the effect of the crack on the stiffness of the structure. The equa- 
tion defining the line spring model is below: 


6 c) 

P 11 P 12 


. P 21 P 22 . 


where « c = displacement due to crack 

e = rotation due to crack 
c 

N = axial force in crack area 
M = bending moment in crack area 

P = crack flexibility matrix. 



(3.2-11) 


The matrix [P] is determined from stress intensity factor calibrations of a 
plane-strain single edge notched specimen using energy-compliance relations. 
The matrix [P] is inverted and added at the appropriate positions in the glob- 
al stiffness matrix. The elastic stress intensity factor is calculated for 
each crack as follows: 


K = (ira) 1/2 [F ^ + F 2 a g ] 


(3.2-12) 
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where K = stress intensity factor 
a = crack depth 
= axial stress 
ag = bending stress 

F 1> F 2 = functions in terms of crack depth and beam thickness. 

The functions Fj and determine the effect of a crack on the stiffness 
of the structure. These variables are functions of the crack depth and beam 
thickness. In addition, these variables are calculated internally in the 
Mechanics of Materials Model (MOMM) and are not required as input. 

The MOMM computer code was designed to be utilized as an initial analysis tool 
for hot section components. Beam elements can be used to create simple three- 
dimensional finite element models that approximate the axial, bending and tor- 
sional stiffnesses of the components being analyzed. The rectangular beam sec- 
tion used in any particular model has stiffnesses that are a function of the 
dimensions of the cross section of the beam. The axial stiffness is dependent 
on the cross sectional area; the bending stiffnesses are related to the mo- 
ments of inertia about the appropriate axes; and the torsional stiffness is 
derived from the polar moment of inertia. Application of simple mechanics of 
materials calculations and engineering judgement are needed to ensure a beam 
design that will produce accurate results. 

3.2.2 Program Validation/Verification 

Some of the test cases (i.e., TEST1 - TESTS) which have been executed to val- 
idate MOMM computer code are summarized below. Each of these cases test vari- 
ous segments of the theory and computer code. 

TEST1 - Cantilever Beam With Axial Load 


A cantilever beam is loaded with a single static compressive loading incre- 
ment. The beam (Figure 3.2-2) is made up of one member, with all degrees of 
freedom constrained at one end and all but two constrained at the end where 
the load is applied. The simple material model is used, and the loading causes 
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only elastic displacements. The lowest frequency and buckling factor are ob- 
tained. The displacements, strains and stresses are found to be: 

u x = P/K = -10' 4 

El = Uj/L = -10" 5 
Oj = Eej = -100 

The resulting lowest frequency and buckling factor are: 

W \/^T 

'lowest ■ jf =Y» - 22.5 

Zir 

Ku 9 

cr = KG’" = 833 ‘ 3 

Agreement between these computed values and independent closed- form solutions 
is exact. 


100 



EA-10 7 

Ku 1 -EA/L - 10 6 

Ku 2 « 12EI/L 3 ■ 10 4 


f 



Figure 3.2-2 Schematic of TEST1 Beam 


36 


TEST2 - Simply Supported, Centrally-Loaded Square Plate 


A quarter of the square plate is modeled using symmetry boundary conditions 
(Figure 3.2-3). Four outside beams and four interior diagonal beams are used, 
with dimensions of the beams chosen so as to reproduce the stiffness and mass 
of the plate. One static loading increment is used with the simple constitu- 
tive model in the elastic range. The nonconstrained degrees of freedom are 
shown . 

The theoretical central displacement is: 

u 2 = .01160 Pa 2 /0 
u 2 = -3.2428 x 10" 6 

The result from the M0W1 computer run is: 


u 2 = -3.4712 x 10" 6 
u 3 



Figure 3.2-3 Square Plate Approximation Centrally Loaded 
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TEST3 - Beam With Axial Enforced Displacement (Static) 


A static analysis (Figure 3.2-4) Is performed using Walker's viscoplastic ma- 
terial model with twelve loading Increments. The properties of Hastell oy X at 

a temperature of 871 “C (1600°F) are used, and the tip displacement is enforced 

-3 -1 

at a strain rate of 3.9 x 10 sec . The computer program reproduces the 
experimental results. A plot of the stress- strain curve obtained from the out- 
put is shown in Figure 3.2-5. 



Figure 3.2-4 Schematic of TEST3 



STRAIN (%) 

Figure 3.2-5 Stress-Strain Response for TEST3 
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ic.3i4 and TESTS - Beam With Axial Enforced Displacements (Transient) 


Both test cases contain a beam fixed at both ends with a node in the middle of 

the beam (Figure 3.2-6). One end is displaced so that the strain rate equals 
3 1 

3.9 x 10 sec" . A transient analysis is performed, with TEST 4 containing 
Walker's viscoplastic material model and TEST5 containing the state-of-the-art 
material model. The viscoplastic material model uses the properties of Hastel- 
loy X at a temperature of 87l’C (1600 # F). Figure 3.2-7 shows the displacement 
at the enforced displacement node, as well as the displacement at the center 
node versus time for each model. The results agree exactly with those obtained 
using a simple Euler integration. 


.3ZAU 

ui At 


t 


.19748 

.00779 


7777 

Figure 3.2-6 Schematic of TEST4 and TEST5 



Figure 3.2-7 Displacement History for TEST4 and TEST5 
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TEST6 - Beam With Crack 


Cracks can be added to the structure at any point between beams and on any of 
the four faces of the cross section. The rules governing the input require 
that at the crack location, two grid points (GP2 and GP3) must be defined with 
identical coordinates. Since the crack is assumed to act at a nodal location, 
there is no length associated with the crack segment; therefore, the user must 
define a vector (v x ) in the direction of the beam length (the element x-axis 
direction). The local y-axis must also be entered (v^) to define the element 
coordinate system. Other data input parameters include geometry set indicator, 
material set indicator, crack depth (D), and a parameter that indicates on 
which face of the cross section the crack exists. Multipoint constraint equa- 
tions must be used to set all the degrees of freedom of the two grid points 
defining the crack (GP2 and GP3) equal to each other except for the axial and 
rotational degrees of freedom. The change in the axial and rotational degrees 
of freedom between the two grid points (GP2 and GP3) is solved for by using 
the crack stiffness and axial force and bending moment at the crack due to the 
loading on the structure (refer to Figure 3.2-8). 



WHERE CRACK 
EXISTS 


Y 



figure 3.2-8 Schematic of Crack Simulation in Beam 
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A beam with a crack one-fifth through its thickness and subjected to a tensile 
load is analyzed. The beam is modeled with two elements. Note that the two 
grid points at the crack, node 2 and node 3, are assigned the same coordi- 
nates. Four multipoint constraint equations are input to set all but the axial 
and bending degree of freedoms equal at node 2 and node 3. 


i: 


EA-10 7 


S 


* 


5 



CRACK 


Figure 3.2-9 


i 


7 / 77777 /- 

Schematic of TEST6 Beam 


The simple material model was used, and the load is in the elastic range. The 
default value of two integration points in each direction of each beam is used. 

The results show that the axial deflection at node 4, which would be 1.0 x 

-5 -5 

10 without a crack, is equal to 1.0204 x 10 . The stress intensity fac- 

tor was found to be 10.931. Using an analytical expression for the stress in- 
tensity factor of a single-edged notch gives a value of 10.860. 


The table below 

shows the comparison between the 

MOMM computer code results 

and analytical 

results for the stress intensity 

factor for different crack 

depths. 


Kl 

K 1 

Crack Depth 

(MOMM) 

(Analytical ) 

.10 

6.442 

6.726 

.15 

8.738 

8.855 

.20 

10.930 

10.860 

.25 

13.374 

13.382 

.30 

16.220 

16.310 

.35 

19.623 

19.818 

.40 

23.764 

23.989 

.45 

28.916 

29.249 

.50 

35.449 

35.845 
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For this example, the value of the K 1 calculated In the MOMM computer code 
agreed with the analytical result except when the crack depth was less than 
0.10. In general, the user should be aware that the accuracy of the calcula- 
tion of K 1 will decrease when the crack depth Is less than one- tenth the 
thickness of the beam. 

TEST7 - Beam With Transverse Hole 


A beam with a transverse hole is loaded with a moment so that the bending is 
in the plane of the hole. The model consists of three beam elements end to 
end, with the middle beam containing the hole. The beam cross section is a 2x2 
square, and the rectangular transverse hole is lxl square. The hole is input 
by setting the transverse hole thickness in the z-di recti on equal to one and 
the middle beam length equal to one. In order to capture the stress in the re- 
gion of the hole, seven integration points in each direction are used in the 
middle beam. The simple material model was used, and the loading is elastic. 



Figure 3.2-10 Schematic of TEST7 Beam 

The displacements, strains and stresses are solved for, taking into account 
the reduced stiffness of the beam with the transverse hole. Interpolating 
using the stresses obtained at the integration points, the stress at the edge 
of the hole is found to be 171.4. The stress at the edge of the beam in the 
section of the hole is equal to 342.9. Stress concentration factors are 
printed out in the output for a circular hole with a diameter equal to the 
thickness of the rectangular hole. The in-plane bending stress concentration 
factor is equal to 2.0, and the stress at the edge of the hole should be mul- 
tiplied by 2.0, i.e., the stress at the edge of the hole is equal to 342.9. 
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List of Symbols 
Referenced Within Section 3.2 


Symbol 

Description 

Page 

1 a| 

Acceleration vector 

30, 31 

Iaf( 

Incremental force vector 

30, 31 

M 

Velocity vector 

30, 31 

CM] 

Mass matrix 

30, 31 

DT 

Time step 

30, 31 

X 

Eigenvalue 

31 

[I] 

Identity matrix 

31 

m 

Eigenvector 

31 
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3.2.4 Problems That Can Be Solved Using MOMM 

The following table presents problems that can be solved using MOMM. 

Material Behavior 


Slmpl Ified 

Walker* s 
Viscoplastic 

State of 
the Art 

Static Analysis 

X 

X 

X 

Transient Analysis 

X 

X 

X 

Lowest or All Frequencies 

X 

X 

X 

Lowest or All Mode Shapes 

X 

X 

X 

Buckling Factor 

X 

X 

X 

Buckling Mode Shape 

X 

X 

X 

Initial or Tangent Stiffness for 

X 

X 

X 

Frequencies and Buckling 
Applied Forces and Moments 

X 

X 

X 

Enforced Displacements 

X 

X 

X 

Thermal Loads 

X 

X 

X 

Cracks 

X 

X 

X 

Holes 

X 

X 

X 

Multipoint Constraints 

X 

X 

X 

Note: Buckling calculations cannot 

be performed in 

a transient problem. 
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SPECIAL FINITE ELEMENT MODEL 


O O 
^ 


3.3.1 Summary 

The finite element code MHOST (MARC Hot Section Technology) has been further 
developed based on the mixed iterative solution technology whose concepts and 
basic hypotheses were defined in the first year effort. These concepts have 
been extended in the current phase to incorporate the effect of multiple em- 
bedded singularities in generic modeling regions. Specifically, a local mesh 
refinement technology has been generated based on the mixed element concept; 
the approach involves a supplemental iteration in conjunction with the intro- 
duction of a higher order polynomial representation for spatial discretization. 

3.3.2 Introduction 

The mixed finite element approximation and its associated iterative solution 
algorithms have been developed for three-dimensional inelastic analyses with 
particular application to turbine engine hot section components. The numerical 
algorithms have been further improved in areas involving accuracy of solution 
and efficiency of computation. 

The enhancement of solution capability has been sought in order to be able to 
deal effectively with problems involving multiple embedded singularities in 
generic modeling regions. The concept of subelement iteration has been derived 
and tested for the present purposes and its numerical performance is shown to 
be superior to that of the conventional finite element method. 

The program development effort includes extensive testing of the capabilities 
built into the MHOST program as well as further enhancements to control the 
iterative procedure in a precise manner. An interface file is generated which 
can be handled by most commercially available finite element postprocessing 
packages. 
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In this section, a number of val Idatl on/veri fl cation problems are Included to 
demonstrate the capability and the performance of the algorithms built Into 
the MHOST program package. 

3.3.3 Literature Survey 

Reports, books and journal articles related to 'nonstandard' finite element 
methods have been surveyed with the main objective being the search for useful 
numerical technology in the framework of the mixed iterative solution approach. 

In the past, 'nonstandard' finite element methods, such as the mixed and hy- 
brid methods, have not been exploited in a systematic manner. The few excep- 
tions generally involve hybrid elements for plate and shell analysis which 
generate element 'stiffness' equations or mixed elements for incompressible 
problems based on the Lagrangian functional used by Herrmann (Reference 1). 

As documented in the First Annual Status Report (Reference 2) and in various 
research papers (References 3 through 6), iterative algorithms for generating 
continuous stress fields such as those developed by Cantin, Loubignac and 
Touzot (Reference 7) can be identified with the mixed finite element method. 
It is important to note that in constructing the iterative methods, use of the 
Hu-Washizu variational principle is crucial to setting up practical useful 
algorithms. In a recent paper (Reference 8), a similar observation was made 
and a finite element method was constructed from the Hu-Washizu principle with 
application to finite deformation plasticity. 

Except for a few early papers such as those mentioned above and literature 
appearing after the present development effort was initiated, no directly 
relevant work has been published describing constructive iterative solutions 
for the mixed finite element equations derived from the Hu-Washizu form. There 
are, however, a number of papers indirectly relevant to and somewhat useful 
for further development of the solution strategy employed in the MHOST code. 
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The algorithm is viewed as composed of three steps. First, the linearized mo- 
mentum equation is solved in terms of the displacement vector for precondi- 
tioning purposes. Second, the postprocessing algorithm is entered to generate 
the strain field based on the mixed interpolation which is used for the inte- 
gration of the constitutive equation. Third, the equilibrium is iterated to 
satisfy the nonlinear virtual work equation with respect to the stress field 
interpolated in a mixed manner. 

It should be noted that the quality of the displacement field generated by the 
preconditioner contributes significantly to the overall quality of the solu- 
tion as well as the convergence characteristics. The use of equivalent stiff- 
ness in the element discontinuous strain mixed forms is a possible numerical 
strategy for these purposes (References 9 and 10). In particular, recent ap- 
plications to plates and shells involving lower order element technology indi- 
cate possible improvements for the preconditioning operations (Reference 11), 
with further efficiency gained by using lower order quadratures (Reference 12). 

The strain recovery algorithm based on the mixed interpolation is found to be 
virtually identical to the classical methodology based on the consistent con- 
jugate stress distribution studied in References 13 and 14. Recent papers 
(References 15 and 16) present a systematic method to construct and analyze 
postprocessing algorithms. It is claimed that there are postprocessing proce- 
dures which provide accuracy for various quantities of the same order as that 
provided by the energy error estimate in finite element displacement algo- 
rithms. This statement agrees quite well with experimental observations by Owa 
(Reference 17) and Nakazawa, Owa and Zienkiewicz (Reference 18). Also, the 
super convergent results, in terms of stress and strain as well as displace- 
ment, often observed in the first year of this project may be due to the 
higher order convergence rate of the postprocessing algorithm used for the 
strain recovery computation. No literature is available on the effect of nu- 
merical quadrature in the postprocessing algorithms except a recent technical 
note by Simo and Hughes (Reference 19) on assumed strain involving the so- 
called B type algorithms for plates and shells. 
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Integration algorithms for nonlinear constitutive equations, in particular 
rate independent plasticity, are well-established as indicated in the Task I 
literature survey (Reference 20); and one of the most reliable algorithms, 
based on the radial return concept (Reference 21) has been implemented. No 
systematic effort has yet been reported on algorithm development for advanced 
viscoplastic constitutive models including temperature effects. Further liter- 
ature search and original investigation is required in this field. 

Methods of equilibrium iteration have been investigated in recent years and a 
number of useful papers and reports have appeared, as mentioned in the pre- 
vious literature survey report. A survey and series of experiments (Reference 
22) provide a useful collection of algorithms and numerical results with ap- 
plication to nonlinear plates. Basic algorithms for the Newton-Raphson and 
Modified- and quasi -Newton methods are compared in the displacement method 
framework. As demonstrated later in this section, algorithms discussed in a 

classical report by Matthies and Strang (Reference 23) are usable even in the 
context of the mixed iterative method. 

The algorithms and convergence arguments directly applicable to the present 

framework are only available from literature on augmented Lagrangian methods 
for the quadratic minimization problem with linear equality constraints such 
as incompressibility and the Dirichlet boundary condition. Possible improve- 
ment of the iterative procedure is indicated in Reference 24, and the numer- 
ical test examples studied in References 25 and 5 show that significant im- 

provement is obtained by using such algorithms for the analysis of Stokes' 
flow. Similar mathematical discussions and algorithms are also found in Ref- 
erences 26 and 27. The original idea of the mixed iterative solution is, how- 
ever, found in a historical work by Arron, Hurwicvz and Uzawa (Reference 28) 
and the class of iterative algorithms for mixed problems is referred to as the 
Uzawa method . 

In the context of linear elastic finite element analysis, the use of mixed 

approximations and equilibrium iterations has appeared in an ad hoc fashion 
repeatedly in the literature other than the work by Cantin, et al (Reference 
7). For nonconforming plate bending elements, Crisfield (Reference 29) has 
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proposed an iterative algorithm to improve numerical accuracy. In the paper, 
the similarity of the concept (based on a modified Hellinger-Reissner princi- 
ple) to the initial strain method iteration is pointed out. Also, an important 
observation by Crisfield is that all significant changes occur in the first 
iteration of the mixed process. 

The implementation of mixed finite elements without an iterative solution is a 
topic that was extensively studied several decades ago, mainly with applica- 
tion to linear problems in mechanics as discussed in a survey by Zienkiewicz 
(Reference 30). These experiences have been used here to avoid possible numer- 
ical difficulties particular to this methodology. This is an important issue 
to be investigated. The characteristics of the particular algorithm used as 
the postprocessor needs to be understood in the mixed method framework, so as 
to fulfill the necessary conditions for stability and convergence. Regardless 
of the solution algorithm, the mixed method needs to satisfy the stability 
condition referred to as the Babuska-Brezzi condition in some sense (Refer- 
ences 31 and 32). However, when equal order interpolations of displacement and 
stress are used, possible violation of this condition is indicated by Oden 
(Reference 33). Resultant oscillations of the numerical solution are indicated 
therein under special circumstances. 

As experienced in mixed/penalty finite element computations for incompressible 
problems, however, the implication of the Babuska-Brezzi condition is not 
quite clear. For instance, as discussed in Reference 34, stability can be 
achieved by using a class of unstable elements which violates the condition ja 
priori but which produces stable results by incorporating a postprocessing 
algorithm which satisfies the necessary condition for the stability. 

Modern development of mixed finite element methods mainly involves construct- 
ing the displacement stiffness matrix from the element discontinuous approxi- 
mation for additional variables (Reference 35). The derivation of this class 
of methods is based on the equivalence theorem stated in Reference 9. An im- 
portant development along this line is a generalization of equivalence theorem 
proposed in Reference 36, indicating that for all numerically integrated dis- 
placement finite element methods, there exists an equivalent class of mixed 
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methods based on the Hu-Washlzu principle. This observation provides some In- 
sights Into the Iterative process with particular application to Inelastic 
problems. As mentioned earlier, these developments of 'reducible' mixed forms 
at the element level are useful mainly to construct the preconditioning system 
of equations as well as to fine tune the equilibrium iterations. 

For the numerical modeling of singularities embedded in a structure, papers on 
fracture mechanics have been surveyed. For displacement finite element tech- 
nology, the report by Fine (Reference 37) covers the technology to date. Ex- 
cept for a few efforts such as reported in Reference 38, not many papers and 
research reports are available discussing mixed and hybrid finite elements 
with application to fracture mechanics, in particular, nonlinear material pro- 
blems. A paper by Babuska and Miller (Reference 16) on postprocessing to cal- 
culate the stress intensity factor was found useful for constructing and 
validating the numerical algorithm implemented into the MHOST program. The 
general strategy for numerical postprocessing is extended to deal with pro- 
blems with singularities. 

A series of papers (References 15 and 16) presents possible utilization of 
postprocessing technology and adaptive mesh refinement. The concepts are 
closely related to mixed iterative solution algorithms in conjunction with the 
subelement calculation as discussed later in this section. Literature on adap- 
tivity and a posteriori error estimate (Reference 39) has been found useful in 
this line of development. 

A series of notes by Axelsson (References 40 and 41) indicate the possible 
utilization of successive relaxation techniques for the solution of finite 
element equations, in particular, the mixed system of equations. However, no 
evidence that such algorithms can be used for nonlinear solutions is provided 
in those publications. 
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3.3.4 Formulation Strategy and Development 


3. 3. 4.1 The Global Solution Strategy 

The nonlinear problem involving the inelastic response of a material body in 
an open, bounded domain fl with sufficiently smooth boundary an is presented 
in this section and an augmented form of the generalized nonlinear Hu-Washizu 
variational principle is derived in the infinitesimal deformation setting. A 
generic solution algorithm is constructed which is valid for quasi -static and 
dynamic- transient analyses. 

The procedure discussed here represents a generalization of the algorithm de- 
veloped in Task I of this project (Reference 2), and the resulting scheme can 
be combined with various solution algorithms and time integration operators 
other than the conventional Newton-Raphson and Newmark-/3 methods. 

The equilibrium equation is given by: 

-<j.. . = p(f. - a-) (3.3-1) 

10,1 J J 

where p is the material density, assumed constant, with £, £ and £ being the 
stress tensor, the acceleration vector and the body force vector respectively. 

We assume a rate constitutive equation: 


a.. = (D t ) L, (3.3-2) 

1J 1 ijkl Kl 

with Dy being the tangent material modulus, and £ and £ the stress and 
strain rates respectively. The strain tensor component is given by: 


:..=4(u- . + u . .)+£?. 
ij 7 i,J J,i ij 


(3.3-3) 


with £° being the initial strain due to thermal expansion and creep effects. 
The usual equality constraints are imposed on the boundary such that: 
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and 



on an 


(k) 

1 


(3.3-4) 


t k = °kl n l = \ 


on an 


(k) 

2 


(3.3-5) 


The weak variational form associated with the above problem statement is: 


(o 1j" u 1j> ’ (p<f j • a j>- u j> * < V V s S ! ^ k) (3.3-6) 

(*«•(*« - (D T ) iJkl 5 kl})= 0 < 3 - 3 - 7 ) 

and 

("ij’ { e 1j ” 7 ( U 1,J * u j.l)}) ’ (°1j’ e 1j ) (3 ' 3 - 8) 

o 

where (...) denotes the usual L inner product over the domain n and <.,.> 
is the integral defined on the specific boundary, with * indicating arbitrary 
variations. 


Elimination of stress and strain from the above system of variational equa- 
tions results in the virtual work equation in terms of displacement: 

a(u, u*) - (f, u*) = 0 (3.3-9) 

with a(.,.) being the usual energy product. The essential boundary condition, 
equation (3.3-4), can be incorporated by virtue of the penalty approach: 

a(u, u*) - (f, u*) + e" 1 < u k - u k , u* > 3n[ k ^ = 0 (3.3-10) 
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or equivalently: 


★ 1 * /b\ 

{o ir u ij ) + e < u k - V u k > 

= (p(fj - a^), ul) + < t k , u k > an!^ 


(3.3-11) 


It can be shown that the simultaneous weak variational statements, equations 
(3.3-6) through (3.3-8), may be derived directly from the Hu-Washizu princi- 
ple. This result implies that the boundary conditions, equations (3.3-4) and 
(3.3-5), enter into the system of equations only via the conservation law for 
linear momentum. In this setting, imposition of boundary conditions for the 
stress/strain mixed mode independent variables is unnecessary. If such con- 
ditions are applied, the well-posedness of the problem may be disturbed yield- 
ing an erroneous solution of possibly a rank deficient system of equations. In 
addition, the penalty approach involving the Dirichlet boundary condition does 
not require the space of admissible displacement variations to fulfill the ho- 
mogeneous counterpart of the same condition. 


Using an equal order interpolation function N for all the variables involved 
in the analysis, we have: 


u k = N M u Mk 
a k = N M a Mk = N M u Mk 
e i j = N M e Mij 
a ij = N M °Mi j 

and for the input (initial) quantities: 

f k = N M f Mk 

ij " n M Mi j 


(3.3-12A) 

(3.3-12B) 

(3.3-12C) 

(3.3-12D) 

(3.3-12E) 
(3.3-12F ) 
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resulting In a system of algebraic equations: 


P 0 B 

M M ^ 

0 D -C 

B T -C 0 

p 

Denoting the L inner products over the domain by (...) and the surface in- 
tegral by the entries of matrices in equation (3.3-13) are: 

P = e" 1 <N T , N> 

B = ( VN, N) 

C = (N T , N) 

D = (N T , D t N) 

M = p (N T , N) 

F = (N T , f) + <N T , t> 

and 

Q = (N T . e Q ) 

Note that, for the sake of simplicity, the time integrated form of the rate 
constitutive equation is assumed. The details of the incremental process and 
the stress recovery are discussed later in this section. Elimination of the 
nodal strains and stresses leads to a 'displacement' solution form: 

B(C)' 1 0(C T )- 1 B T + P| u = F - M a (3.3-14) 
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with 


F = F + B(C) T DttV 1 Q, (3.3-15) 

whereas the standard finite element result based on equation (3.3-9) is ap- 
proximated by a somewhat simpler form: 

(K + P) u = F - M a + Q\ (3.3-16) 

with Q* being the nodal force generated by the initial strain in the structure. 

An iterative solution algorithm for the quasi -static counterpart of equation 
(3.3-13) can be constructed as follows: 

(a) Set a vector R = 0; initialize the displacement vector u = 0 

(b) Solve the preconditioning equations to update the displacement such 
that: 

u = u + A" 1 (F - R) (3.3-17) 

(c) Recover the nodal strain: 

e = (C T ) _1 (B T u - G) (3.3-18) 

(d) Integrate the constitutive equation: 

a = £ D t e dt* (3.3-19) 

with t* being the quasi-time associated with deformation history. 
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(e) Evaluate the residual: 

R = B T s (3.3-20) 

(f) If the residual is small enough, then exit; or else repeat from step 
(b). 

As is obvious from the above discussion, the choice of the preconditioner A is 
crucial in obtaining convergence characteristics necessary for a practical im- 
plementation. For instance, if: 

A = B(C)' 1 D(C T )- 1 B T + P, (3.3-21) 

then no iteration is needed. However, sparseness properties of the finite ele- 
ment system matrix could no longer be exploited if the above form were to be 
utilized. As a reasonable compromise, we use the equilibrium equation in- 
volving an augmented displacement stiffness matrix which is set to: 

A = K + P. (3.3-22) 


Other methods of preconditioning have been investigated, but so far no robust 
scheme is known to be applicable for a wide range of solid and structural 
analyses (Reference 42). Except for minor modifications, the present solution 
utilizes the form defined by equation (3.3-22). 

3. 3. 4. 2 Incremental Iterative Solution Algorithms 

The inelastic problem is solved through the deformation history in an incre- 
mental manner. Let the solution for a given state of displacement, strain and 
stress which satisfies the nonlinear algebraic equation (3.3-13) for specified 
load and displacement boundary conditions be identified by u 11 , je n , s 0 , 
and F n , respectively. Then for a given load increment, AF n , the increments 
of displacement, strain and stress denoted by aju 11 , Aje n , and as 0 are de- 
termined in an iterative fashion as described in the previous section. Note 
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that the tangent stiffness Ky is used instead of the total stiffness equa- 
tion in the preconditioning process. To follow a complicated equilibrium path, 
an automatic adjustment procedure which controls the size of the load incre- 
ment in the iterative process is available. The algorithm in outline form is: 


(a) Set the residual vector R = 0, initialize the incremental displace- 
ment vector au 11 = 0. If the first increment, initialize the total 
load factor X. 


(b) Project the displacement by: 


du = 



(3.3-23) 


and if it is the first cycle of iteration, calculate the arc length 
and initialize the incremental load factor (ax) and dx. The arc 
length l is defined by: 


* = 


(du) 


1/2 


and the incremental load factor is given by solving the quadratic 
equation for the iterative load change. Then, carry out the back sub- 
stitution to get: 

du = Kj 1 ( XF - R) (3.3-24) 

(c) Update the displacement vector: 


au = au + du + dX du 


(3.3-25) 
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(d) Find the total load factor update: 

X = X + dX (3.3-26) 

and the incremental load factor: 

AX = AX + dX (3.3-27) 

based on the spherical path formulation (Reference 43). 

(e) form the residual vector in the mixed manner [equations (3.3-18 
through 3.3-20)] and then check the convergence. If convergent, start 
the next increment; or else repeat step (b). 


Figure 3.3-1 presents a flow chart for this algorithm, and further details re- 
garding this process will be included in the MHOST Theoretical Manual. In the 
mixed method, steps which generate the residual force are treated as a package 
of operations. Any iterative method designed to improve convergence character- 
istics can be employed. For instance, the BFGS update procedure (Reference 23) 
is utilized in the following fashion: 


(a) Initialize the residual vector JR and the incremental displacement aij 
such that: 


R = 0; au = 0 


(3.3-28) 


(b) Modify the residual with i being the iteration counter: 


• 'I T f 

5-13-2 V)? 


(3.3-29) 
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{ c ) I ntermedi di spl aCcmcivt update* 


du = 



(3.3-30) 


where K-J-"* is the latest stiffness matrix factorized in the current 

increment. The BFGS iteration counter is set to 1 when the factori- 
zation is performed. 


(d) Complete the displacement update: 

I! 


d “ ■ (3.2 (I * ?j !/> 


du 


(e) Form the new incremental displacement: 
au = au + du 


(3.3-31) 


(3.3-32) 


(f) Form the new residual with respect to the updated incremental dis- 
placement using equations (3.3-18) through (3.3-20). 

(g) Cneck convergence and, if necessary, repeat from step (b) or else, 
exit. 

Note in the above algorithm, vectors v^ and represent the iterative 
changes in the residual and displacement vectors respectively so as to form 
the inverse BFGS update: 

[K Ti ] _1 = (I + w. vT) LK^]" 1 (I + v. w. T ) (3.3-33) 
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Figure 3.3-1 The Arc-Length Method for the Mixed Iterative Solution 

Figure 3.3-2 presents a flow chart for the BF6S procedure. It is possible to 
introduce combined BFGS - arc length algorithms in the mixed iterative solu- 
tion algorithms incorporating a line search technique. The applicability of 
such advanced techniques requires further investigation. 
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figure 3.3-2 The Inverse BFGS Update for the Mixed Iterative Solution 



3. 3. 4. 3 Construction of the Stiffness Equations 


A number of different element based integration schemes are implemented into 
the MHOST numerical solution procedures. To achieve stability and coarse mesh 
accuracy in the finite element solution, not all the terms in the equations 
are integrated at the nodes. Two-point Gauss quadrature is used for evaluation 
of the stiffness coefficients and the internal force vector in each coordinate 
direction. A nodal quadrature {element discontinuous trapezoidal integration) 
is used in the strain recovery operations, and all the stress components are 
evaluated at the global nodes. 

To extract near optimal numerical performance of the elements implemented in 
the MHOST code, a selective reduced integration option is available in which 
the shear strain components are evaluated at the internal element centroid. 

All the integration procedures are internal operations designed to generate an 
accurate finite element solution and are invisible to the user as results are 
reported primarily at the nodes. 

Two major refinements have been incorporated into the MHOST program to ensure 
the generation of a good displacement update for preconditioning purposes. 
These are an improved version of the filtering scheme for selectively reduced 
integration and a modified numerical quadrature for plates and shells which 
avoids possible kinematic mode excitation. 

Coordinate Transformation and Filtering Algorithms 

With application to general two- and three-dimensional elements, a method is 
developed for construction of the element coordinate system and its utiliza- 
tion to filter particular strain components for the selectively integrated 
stiffness equations. A nonstandard notation is used to maintain maximum pos- 
sible generality in the following discussions. Whenever indicial notation is 
applied, the summation convention is assumed for repeated indices unless 
otherwise stated. 
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In this section, we first establish the coordinate systems used in the iso- 
parametric finite elements. In particular, the local (pointwise) definition of 
the orthogonal system associated with the Jacobian matrix is presented based 
on the theory of polar decomposition, and this definition is used later to de- 
vise a family of filtering algorithms for the element strain components. The 
deformation tensors associated with the Jacobian matrix are used as a basis 
for constructing the measure of isoparametric distortion. The filtering scheme 
is constructed from strain tensor components with particular application to 
selectively reduced integration to avoid possible numerical locking and fur- 
ther to improve the accuracy of finite element displacement type solution 
procedures. 

Consider a linear isoparametric element in two dimensions. We denote the glo- 

1 2 

bal coordinates fixed in physical plane by x = (x , x ) and the isopara- 

~~ 12 

metric coordinates in the reference plane by £ = (q , q ). The notation 
and the results which follow are readily extended to general three-dimensional 
elements. 

We denote the Jacobian matrix by J^, which is defined by: 

J = [J k K ] (3.3-33) 

where 

J k K = x k , K = 3 x k /aq K (3.3-34) 

By definition the Jacobian is invertible and for any invertible linear trans- 
formation there exists a decomposition, known as the polar decomposition such 
that: 


J = R U, 


(3.3-35) 
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where R Is the orthogonal rotation tensor and U Is a symmetric, positive- 
definite matrix often referred to as the right stretch tensor . In the present 
setting, the rotation tensor represents the orientation of the physical coor- 
dinate system with respect to the orthogonal basis associated with the iso- 
parametric coordinates in physical space. Hence, its inverse gives the basis 
vectors for the orthogonal coordinates parallel to the element orientation. 

The computational procedure to obtain R“* is to utilize the relation: 

J T J = U T U, (3.3-36) 

as demonstrated in Reference 44 among several others, together with the decom- 
position, equation (3.3-35): 


R = J U" 1 (3.3-37) 

The inverse of the right stretch tensor is obtained from the eigenvalues and 
eigenvectors of the right Cauchy-Green tensor C which is defined as: 

C = J T J. (3.3-38) 

Hence, 

IT 1 = -1/2 = nJ *‘ 1/2 N c (3.3-39) 

C 

with and being the eigenvalues and eigenvectors, respectively, ar- 
ranged in matrix forms. Hence, the equation used in the actual computation is: 

R = J N c T X~ 1/2 N q . (3.3-40) 

The purpose of selective integration is to avoid numerical locking due to the 
overestimation of certain element deformation modes as illustrated in Refer- 
ence 34. This is realized by integrating the contribution of strain energy 
associated with the strain components characterizing such deformation modes at 
the reduced integration points. The strain components to be under- integrated 
are not defined in the global coordinate system but in the coordinate system 
parallel to the element orientation given by the rotation tensor R discussed 
above. 
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Let e-^- be the global Cartesian tensor associated with the strain 

tensor component in the local coordinate system. Denoting the coordinate 
transformation matrix by g such that: 


g = [gjj] = R _1 (3.3-41) 

where the superscript and the subscript denote the local and the global coor- 
dinates, respectively. Then we have: 


S kl) 

ij 


,k i ki 
J i g j ‘ 


e*n* = gv g,- e- = g- 9? g; g; e 


k 9 k 9 
i ^m ^n °mn* 


(3.3-42) 


Introducing an array G referred to as the filtering matrix such that: 

Ic k k 

G^. = g.j gj (no sum on k), 

we simplify the above expression as: 

t < k ‘> = G k 6* « 
ij im jn mn 

Hence, the volumetric strain in the element coordinate system is given in 
terms of the global strain components by: 


(3.3-43) 


(3.3-44) 


(v) r k r £ ki 
e . . = b . b . 6 e__ 

ij im jn mn 


and the deviatoric strain by: 

e<°>. -SL G< (1 - « kl ) 


ij im jn 


'mn 




where 6 is the Kronecker delta. 


Defining the original strain component e by: 

e = B u, 


(3.3-45) 


(3.3-46) 


(3.3-47) 


where B is the usual strain- displacement matrix calculated at the quadrature 
points and u is the nodal displacement vector, the energy functional I(u) is 
obtained in terms of nodal displacement for linear elasticity as an example by, 

I( u ) = y\ e (v) :e (v) dx + e (D) :e (D) dx - Jf * u dx - Jl ’ u ds (3.3-48) 

n fi n as? 
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where : denotes tensor Inner products, X and y are the Lame-Navi er coeffi- 
cients and f and T are the prescribed body force and surface traction respec- 
tively. 

In the above energy functional, either the first term or the second term is 
under- integrated to obtain the necessary effects of the selectively reduced 
i ntegrati on . 

In computations, it is often more convenient to construct the strain-dis- 
placement matrix selectively, so that the above energy functional can be writ- 
ten as: 


I(u) = u T / bI D B- dx u - F T u (3.3-49) 

n o o 

where F is the collection of prescribed load terms which appeared in equation 
(3.3-48). The new strain-displacement matrix B c is the selectively sampled 
matrix equivalent to a certain mixed method under isoparametric distortion. In 
matrix notation, the filter for the element volumetric strain is written as: 

e (v) = G T I G B u (3.3-50) 

and for the element deviatoric strain: 

e (D) = G T (1 - I) G B u (3.3-51) 

at each integration point, where 1 and I are the matrices: 

1 = [1] ; I = [6.j] . (3.3-52) 

Hence, we have: 

B $ = G T I G B + G r (1 - I) G B (3.3-53) 

where either the first or the second term is sampled at the reduced integra- 
tion points and substituted to the array associated with the full integration 
points. Such an operation is trivial only for the linear Lagrangian elements 
with a single reduced integration point. 
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Hourglass Control Algorithm for Plates and Shells 


The kinematic model generated by the reduced integration of transverse shear 
terms in 4-noded bilinear plate and shell elements often induces numerical 

noise in the displacement preconditioning operation. However, excitation of 

these modes, often referred to as the hourglass modes , does not cause signi- 
ficant deterioration in the convergence of the mixed iterative algorithms. To 
eliminate these modes, filtering methods have been developed either to con- 
strain the modes by modifying the stiffness equation ( a priori hourglass con- 
trol) (Reference 45) or to filter out the sprious noise after the nodal dis- 
placement is obtained ( a posteriori hourglass control) (Reference 46). 

A simplified a priori hourglass control algorithm based on an approach similar 
to the scheme proposed in Reference 47 has been implemented into the MHOST 

code. The algorithm takes advantage of the fact that the fully integrated 

stiffness matrix (2x2 integration) does not contain kinematic modes. 

The transverse shear stiffness matrix K is constructed as follows: 

~s 

K $ = e K* 2x2) + (1 - e) K< lxl) (3.3-54) 

where e is a small parameter associated with the aspect ratio of the element. 
Numerical locking is avoided due to the insignificant participation of the 
fully integrated terms. 

The value for e is calculated at the centroid of the element using the formula: 

e = e 0 t c / h (3.3-55) 

where t is the thickness of the element at the centroid and h is the mesh 
c 

size given by: 


h 



- x 4 | + | x 2 



(3.3-56) 
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Typically a value of 0.01 for e Q has been used for the validation and veri- 
fication exercises. 

No additional computational cost is involved in this hourglass control algo- 
rithm because the element stiffness equations are integrated selectively and, 
hence, K^ 2x2 ^ is available without adding a new integration procedure. 

3. 3. 4. 4 Time Integration Algorithm 

The basic temporal discretization procedure is virtually identical to the 
quasi -static solution of the mixed finite element equations. The system of 
ordinary equations: 


M ii = I (u) (3.3-57) 

where I is the nonlinear function of nodal displacements and is discretized in 
time resulting in a recursive form for updating u such that: 

A t u t+At = B* (u*) (3.3-58) 

where A and B are the linearized versions of the time integration operator 
matrices. The final goal here is to use formula (3.3-58) iteratively as the 
preconditioner to find a displacement vector and its time derivatives which 
sati sfy : 


I (u t+at ). (3.3-59) 

The generalized Newmark-/3 method is used to construct the iterative solution 
algorithm. The semi-discrete approximation yields a set of nonlinear ordinary 
differential equations: 


Ma + Cv + / B T a dx = F, 
Q 


(3.3-60) 
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where the third term on the left-hand side is the integral 
strain energy, and for elasticity we can write: 


associated with the 


JV O dx = K u . (3.3-61) 

ft 

Otherwise, the deformation history needs to be integrated for evaluating this 
term. The mass matrix is defined in terms of the nodal basis functions, and 
mass density is designated by the symbol p: 


M = / p N T N dx (3.3-62) 

ft 

One of the possible forms for the damping matrix is to express it in terms of 
the mass and stiffness matrices as follows: 


C = Cj M + c 2 K . (3.3-63) 

Based on the weighted residual argument of Zienkiewicz (Reference 30), the 
generalized form of the Newmark-/3 algorithm is written in a recursive manner 
as: 


u t+ At = u t + Atv t + At^ | (1 _ 2i3 ) a t + z(3 a t + a 1 1 (3.3-64A) 

U t+At = v t + At | (i . y ) ^ + Y a t+At | (3.3-64B) 

where t indicates the current time (known) and At is the current time incre- 
ment. 


The overall equilibrium at the next time level is: 

Ma t+At + Cv t+At + f B T a (u t+At ) dx = F t+At . (3.3-65) 

ft 

Approximating the energy term linearly yields: 

f B T a (U t+At ) dx = K* AU + f B T a (u 1 ) dx (3.3-66) 

ft n 
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where K* Is the tangent stiffness matrix which shall be considered later. 
Equilibrium at the current time level: 


implies that: 


Ma* + Cv* + f B T a (u*) dx = F* 


[b T a (u t+At ) dx = K* au + F t - (Ma* + Cv*) . (3.3-67) 

Substituting equation (3.3-67) into (3.3-65), we have the following incremen- 
tal equilibrium form, where au is the displacement increment associated with 

At: 


Ma t+At + Cv t+At + K* au = aF* + (Ma t + Cv*) . 
From equations (3.3-64A) and (3.3-64B): 

t + At 1 ... 1 w t /I n -t 

V“'P ' ( ^‘ 1 


(3.3-68) 


(3.3-69A) 


v t+At = _i_ au - d - - 2 j|) v* + At [(1 - t) - t (^ -1)] a 4 . (3.3-69B) 

Substituting the above, equations (3.3-69A) and (3.3-69B), yields a linearized 
algebraic system of equations: 


1 v t + 

j8At 


^ M + C + K*J au = aF* + (Ma t + Cv t ) + M | 

(^ - 1) a 1 | + C | (1 - ^) v - At [(1 - y) - y - 1)] a t | 


(3.3-70A) 


which we write simply as: 


A A 

K AU = F 


(3.3-70B) 


Now, we shall recast the above algorithm in the mixed iterative framework. 
First, the conventional Newton-Raphson scheme is constructed in the context of 
the displacement calculation: 


A 1 A A 

A Vl = AU n + V (F - V 


(3.3-71) 
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with R n beina the residual at nth iteration cvcle defined bv: 

II V ' ■ -V - - - — V - 

R n = /B T { s* + as (au 11 ) } dx - |Ma t+At + Cv t+At } . (3.3-72) 

The updated values for acceleration and velocity vectors in equation (3.3-72) 
are calculated from the formulae (3.3-69A) and (3.3-69B) using the latest up- 
date for incremental displacements. The mixed interpolations for stresses are 
used in place of total and incremental stress arrays in the first term of 
equation (3.3-72). The computation is accomplished in exactly the same manner 
as in the quasi-static incremental iterative analysis. 

Introducing a new variable d R : 

d n = K" 1 (F - R n ), (3.3-73) 

the iterative update of the velocity and acceleration is written as: 

a^ At = a^ At + -L d (3.3-74A) 

n+1 n n 

v£ At = v** At + d n (3.3-74B) 

n+1 n ySAt n 

with the starting values: 

(3 - 3 - 74C) 

vJ +At = - (1 - ^) v t + At [(1 - y) - y (^ - 1)] v 1 (3.3-74D) 
calculated at the beginning of each time increment. 

Note that the load vector F is left unchanged throughout the iteration. 
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3. 3. 4. 5 Eigenvalue Extraction for Modal and Buckling Analysis 

The eigenvalue extraction procedure In the MHOST program utilizes a standard 
subspace Iteration method to obtain the natural vibration frequencies and the 
buckling load at a given quasi -static load step. The band matrix solver Is 
used to factorize the global stiffness equations prior to the subspace itera- 
tions. In the subspace, the Jacobi method is used to extract all the eigen- 
values and eigenvectors. 

For a modal analysis, the consistent mass matrix is formed including the rota- 
tional inertia terms for shell elements. The generalized mass is calculated 
and reported in conjunction with the eigenvalues and eigenvectors. 

In a buckling analysis, the initial stress matrix is calculated from the nodal 
stress resulting from the mixed iterative solution for the initial quasi - 
static loading. The displacement stiffness equations are used to represent the 
structural model in the eigenvalue analyses. However, the accurate represen- 
tation of the initial stress terms Improves the overall accuracy of the buck- 
ling load estimates. 

3. 3. 4. 6 Subelement Iterations 

Computational fracture mechanics aspects of the inelastic analysis of turbine 
engine hot section components are discussed in this section. The motivation 
for this endeavor is to seek an economically feasible numerical process with- 
out sacrificing the accuracy of the solution. 

The standard finite element method is to take into account the effect of em- 
bedded singularities by refining the mesh subdivision In the neighborhood of 
such points, or alternatively to Introduce special elements with singular 
functions In this neighborhood. The currently available approaches are often 
prohibitively expensive, especially when the analysis Involves a structure 
with multiple singularities, each of which needs special treatment. 
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The use of nonstandard schemes, in particular the mixed finite element form, 
is accurate and stable, and appears to be far more advantageous when compared 
with conventional schemes in regular nonlinear structural analysis. Hence, 
utilization of nonstandard schemes for problems including singularities is 
worth investigating to see if a major improvement in computational fracture 
mechanics can be realized. 

The main reason for using the mixed finite element method in problems with em- 
bedded singularities is that higher resolution for stresses and strains can be 
obtained without excessive mesh refinement for the displacement variables near 
the singularities. 

In the framework of the displacement finite element method, the only informa- 
tion primarily available from the computation provides the nodal displacements 
from which the strain/stress components are calculated by differentiating the 
shape functions. The stress/strain values are then used to determine fracture 
mechanics related quantities such as the stress intensity factor and the J- 
i ntegral . 

It is clear then that the accuracy of the strain/stress components is a full 
one order less than that of the displacement components for the displacement 
finite element method, even without a singularity. Moreover, the approximation 
of strains/stresses could be extremely unstable near singularities, as well as 
in regions where the deformation is highly concentrated. To obtain accurate 
results for a problem with embedded singularities by the standard displacement 
technology, mesh refinement near singularities is unavoidable in order to be 
able to determine an accurate and stable displacement field, which indeed is 
the only source for generating the strain/stress approximations. For the best 
possible results, monitoring an a posteriori error indicator, together with 
several passes of (adaptive) mesh refinement are advisable in such calcula- 
tions. 
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When the strains or stresses, sometimes both, are Included explicitly In the 
finite element system of equations. In particular Interpolated by (^-con- 
tinuous basis functions, the resulting approximations are not only accurate 
but also stable, having no spurious modes In the strain/stress recovery oper- 
ator. 

For regular problems, the numerical representation of the stress field in par- 
ticular is quantitatively accurate due to the equilibrium condition being ex- 
plicitly satisfied by the approximate stress field itself. Compared with dis- 
placement method results, a full one order improvement is indeed expected of 
the convergence rate in the mixed method for the same mesh subdivision. It 
should be noted that the role of the displacement solution in the mixed form 
is to generate qualitatively the deformation mode from which indirectly the 
deformation gradient is extracted and fed into the stress recovery and equi- 
librium iteration operations. Therefore, the quantitative measure of error in 
the displacement vector contributes relatively less significantly to the over- 
all error in calculating fracture related quantities. 

The approximate solution procedure for a deformable body with embedded singu- 
larities involves first solving the total structural problem without the sin- 
gularities; this step shall be specifically referred to as the global system 
approach. The second step then computes the deformation and the stresses near 
the singularity based on the first step results and is called the local system 
approach. The concept is somewhat similar to substructuring in conventional 
finite element computations. In Figure 3.3-3, the heirarchical structure of 
the mesh subdivision and subelement representations is illustrated. 

Assume that P in Figure 3.3-3(a) is a singularity embedded in the structure 
the size of which is small compared with the scale of the whole structure. 
Typically, the conventional finite element representation needs a quite fine 
mesh for such a situation. Figure 3.3-3(b) is an example of a reasonable mesh 
subdivision to represent the global behavior of the structure. In this model, 
the singularity is still embedded in an element, but not explicitly considered 
in the numerical model. This model is referred to as the global finite element 
mesh. 
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The subelement mesh subdivision Is shown In Figure 3.3-3(c) where the element 
with an embedded singularity is subdivided further down to the scale of the 
singularity and its effect is now evaluated explicitly. 



A SINGULARITY EMBEDDED IN 
A STRUCTURE 



THE SINGULARITY EMBEDDED IN 
AN ELEMENT 



THE SINGULARITY REPRESENTED 
BY THE LOCAL SUBELEMENT 
MESH 


Figure 3.3-3 Local -Global Representation of an Embedded Singularity 
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Topologically, the global mesh Is constructed such that the minimum element 
size Is much larger than the scale of the singularities so that no singularity 
is identified by the global system. The local system Is constructed In a 
single global element and is coupled with the global system through nodes com- 
mon between the single global element and the local elements within it. 

A major difference of the present approach from substructuring is that the 
local refinement takes place a posteriori and the information extracted from 
the locally subelemented element is brought back to the global mesh via the 
residual load correction. This implies that no special algebraic treatment is 
required even when material nonlinearity occurs in the subelement region. 

The concept of interacting two geometrical models implies an iterative algo- 
rithm as follows: 

(a) Solve the global system of finite element equations for precondi- 
tioning purposes, equation (3.3-17). 

(b) Use the above results together with information available for the 
singularity to set up the local finite element equations which have, 
in principle, the same form as equations (3.3-17) through (3.3-19). 

(c) Form the internal residual load vector based upon the result of step 
(b) and modify the right-hand side vector of the global system. 

(d) If the residual load vector is small enough, then exit; otherwise re- 
peat from step (a). 

As a result, overall equilibrium is achieved with respect to the presence of 
the local system. 

This solution strategy represents a mixed version of the adaptive process 
which has been investigated in recent years and which has yielded results that 
are encouraging for elastic stress analyses with singularities. The computa- 
tional strategy we propose here is to define the local mesh and the refined 
interpolation altogether in the region near the singularity in a manner some- 
what similar to that of combined h-p refinement. 
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As a result, we expect to obtain an accurate stress-strain solution in the 
global approximation subspace enriched by the local approximations. 

The final solution is found in the global -local system with the global system 
being used only for preconditioning purposes. This means that factorization of 
the approximate global stiffness matrix obtained by the conventional displace- 
ment method may be sufficient to kick off the present solution procedure. Fur- 
ther details will be discussed in following sections. 

Consider as another practical situation a single displacement element with an 
embedded hole. The isoparametric element in the physical plane shown in Figure 
3.3-4(A) is mapped onto the reference plane as shown in Figure 3.3-4(B). 

Some computational aspects need to be considered here for application of the 
present method to practical problems of turbine engine hot section components. 

From a programming point of view, a parametric representati on of an elliptic 
hole by its size, orientation and location in a displacement element is con- 
venient because this allows the code to generate directly subelement mesh data 
whenever such data are needed. The data structure is simple as no permanent 
storage allocation is required for the subelement mesh. In Figure 3.3-4(B), an 
example of such a ready-made subelement mesh is presented for a circular hole 
located at the center of displacement element mesh. 

On the other hand, it may be user-friendly and perhaps conceptually more gen- 
eral to define explicitly the subelement mesh as an additional data set. The 
data structure then needs to be reviewed so as to allocate additional memory 
for subelement mesh storage. 

As experienced in the h-version of adaptive mesh refinement, it is anticipated 
that the subelement mesh approach will be used recursively in a fashion simi- 
lar to that of the multi-level substructure technique. A well -organized data 
storage scheme needs to be developed, therefore, in order to realize a fully 
flexible implementation of the proposed method. 
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figure 3.3-4 Multiple Isoparametric Mapping for an Element with an Embedded 



In the following discussions, we assume (for the purpose of simplicity) that 
the subelement mesh Is defined in the isoparametric element coordinate system 
rather than in the physical space in which the global mesh is defined. 

The stresses and strains are represented at the nodes of subelements defined 
in the first reference plane. Figure 3.3-4(B). We introduce variables to char- 
acterize the internal deformation of the stress and strain subelements so as 
to obtain these quantities uniquely in an accurate manner. As mentioned in the 
first section, some further mathematical investigation is required to come up 
with an optimal combination of subelement mesh definition, functions to repre- 
sent the deformation of each subelement as well as the stress-strain inter- 
polations. 

A major difficulty encountered in multidimensional computation involves inte- 
grating the coupling matrix B over the subelement. Figure 3. 3-4(0. Once the 
size of the hole and its location are identified, the values of displacement 
interpolation functions are obtainable at every nodal point of the subelement 
mesh. Denoting these quantities by N K , we introduce the local displacement- 
type interpolation of N in these elements by: 

N = N k S N k (3.3-75) 

and the coupling matrix is integrated approximately by: 

f aN? c 

“ = [ s¥rbT H L df d " \ ] - (3 - 3 - 761 

where the superscript S denotes the functions defined at the subelement level. 
The calculation is carried out on each local subelement 0 S and summed over 
the global displacement element. This general integration procedure plays a 
central role in coupling the global and local meshes by representing the resi- 
dual load vector at the global displacement node locations yet reflecting the 
information at the local element level. For consistent definition of the oper- 
ator matrix B, the transformation of the shape function from global to local 
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needs to be carried out as defined by equation (3.3-75) on the Isoparametric 
element space. Figure 3.3-4(B). This must also be properly coordinated with 
the fact that each subelement In Figure 3.3-4(B) Is Individually mapped Into 
another Isoparametric space. Figure 3.3-4(C). 

It is obvious that the location of nodal points Is needed in the element coor- 
dinate system. When the locations of subelement nodal points are specified in 
the physical plane, then the inverse of the nonlinear Isoparametric mapping 
needs to be calculated in order to find these points in the element coordinate 
system. 

We outline the solution strategy in a general format. First, as a precondi- 
tioner, we introduce the conventional stiffness equation with respect to the 
global displacement degrees of freedom u as follows: 

K u = F . (3.3-77) 

This is derived directly from the virtual work principle in terms of displace- 
ment. A modified recursive form of the mixed finite element equations is writ- 
ten again in terms of the global unknown variables: 



Setting ju^ 0 ^ =0 and =0 with the superscript being used for the 

iteration counter, we solve the system of displacement equations: 

u (n+l) = u (n) + K -1 (F _ Ba (n) )< 


(3.3-79) 


Then the stresses and strains are updated in the subelement region by solving 
implicitly: 



with u g being the enriched displacement in the subelement region. To obtain 
this quantity, the factorization of a displacement type stiffness equation is 
required in this region. 

It should be noted that the discussion presented here describes in general 
terms the solution strategy used. The approach actually implemented into the 
code involves hierarchical displacement approximations for realizing the in- 
teraction between the global and local meshes. 

Equation (3.3-80) is improved by increasing the information contained on the 
right-hand side. The use of additional terms characterizing the deformation 
near the singularity needs to be considered. The simplest approach, for ex- 
ample, is to take the displacement as a dummy variable at stress- strain nodal 
points, with the displacement interpolated by conventional polynomial basis 
functions in a manner similar to that of the substructuring technique. 

When such variables are introduced, the solution procedure becomes very close 
to what is called the multi -grid method in the finite difference context. The 
overall solution flow can then be given as follows: 

(a) Solve for the global displacement by equation (3.3-79). 

(b) Recover the strain at nodes in the global finite element mesh. 

(c) Evaluate the stress at nodes in the global finite element mesh. 
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(d) Enter the inner loop and calculate the displacement quantity (used 
only as a dummy variable) in the subelement mesh. 

(e) Evaluate the stress and strain by equation (3.3-80). 

(f) Check the convergence in terms of residual in the subelement mesh. If 
not converged, repeat from step (d). 

(g) Evaluate the global residual including the stresses in the subele- 
ment. If convergent, start a new increment; otherwise, repeat from 
step (a). 

Through the inner loop, the global finite element mesh interacts with the sub- 
element and an accurate, high- resolution stress field is obtainable without 
increasing the size of the global stiffness equations. 

It is emphasized that additional coding is needed to incorporate steps (d) to 
(f), which can be performed by adding an element routine to handle the em- 
bedded singularities and by allocating the core storage associated with the 
arrays used for the inner iteration. As the size of the algebraic equations 
used for representing the embedded singularity is limited, a performance test 
is needed to decide whether to employ either a nested loop for within-element 
solution or factorization of the same matrix without iteration. 

3.3.5 Computer Program Development 

3.3.5. 1 Solution Capabilities 

Version 2.0 of the MHOST program, which is the development version for Task II 
efforts, currently supports many options and a limited number of linear and 
quadratic finite elements, all of which are operational for both the mixed 
iterative approach described in Section 2 and the conventional displacement 
method. The MHOST analysis capabilities are summarized in Tables 3.3-1 and 
3. 3-1 I. Additional information is provided below through descriptions of the 
control parameter keywords. 
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Table 3.3-1 

MHOST Analysis Capability: Element Definition Options 


Element 

Definition 

Options 


Axl- Three- Three- 

Plane Plane symmetric Dimensional Dimensional 

Beam Stress Strain Solid Solid shell 


Linear** 

Elasticity x x x 


Simplified 

Plasticity x 

*1 

Elasto- 

Plastlclty x 


Unified 

Creep-Plasticity 

*2 

Thermal 

Strain 


. *2 
Creep 

Strain 


NOTES: 

*1 Applicable to Isotropic and anisotropic (user subroutines) 
materials. 


*2 Not applicable to unified creep-plasticity In which these 
quantities are the Integrated part of the model. 


x*x x x 

xxx x x 

xxx x x 


xxx x 

xxx x 


Table 3.3-II 

MHOST Analysis Capability: Analysis Module Options 


Analysis 

Module 

Option 

Beam 

Plane 

Stress 

Plane 

Strain 

Axl- 

symnetrlc 

Solid 

Three- 

Dimensional 

Solid 

Three- 

Dimensional 

Shell 

Quasi- 

static 

Analysis 

X 

X 

X 

X 

X 

X 

Buckling 

Analysis 

X 

X 

X 

X 

X 

X 

Nodal 

Analysis 

X 

X 

X 

X 

X 

X 

Transient 

Dynamics 

X 

X 

X 

X 

X 

X 
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♦ANISOTROPY 


Elastic material anisotropy may be included in an analysis by adding the user 
subroutine UHOOK. The anisotropic plastic response of a material is described 
by the user subroutine ANPLAS as documented in the MHOST User's Manual. 

♦BFGS 

The inverse BFGS update procedure is invoked by flagging this optional parame- 
ter. The default interative algorithm is the conventional Newton-Raphson 
scheme. 

♦BOUNDARY 

Nodal displacement constraints are imposed using a penalty approach as dis- 
cussed in Section 3. 3. 4. 2. 

♦BUCKLE 

The initial stress stiffness matrix is formed and an eigenvalue extraction is 
performed to obtain buckling modes. This option can be invoked at an arbitrary 
step of the incremental nonlinear solution process in order to detect the 
change in the buckling load due to inelastic response of the structure. 

♦CONSTITUTIVE 

Three different constitutive formulations are included in the code for de- 
scribing material behavior. They involve: (1) secant elasticity (simplified 
plasticity) in which the material tangent is generated for use with Newton- 
Raphson type iterative algorithms; (2) von Mises plasticity with the asso- 
ciated flow rule treated by using the radial return algorithm; and (3) the 
nonlinear viscoplastic model developed by Walker, in which an initial stress 
iteration using the elastic stiffness is utilized. A linear elasticity option 
can be flagged for experimental purposes, and the default is the conventional 
von Mises plasticity model. Anisotropic plasticity is invoked by updating the 
user subroutine ANPLAS. The dummy subroutine supplied with the MHOST code is 

m 

always required to produce correct results for isotropic cases. 
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♦CREEP 


Creep effects are taken into account by integrating the time history in an ex- 
plicit manner. An optional self-adaptive time step size control algorithm is 
available. 

*D I SPLACEMENTMETHOD 

This option invokes the conventional displacement model in which the residual 
load is evaluated directly at the integration points. For linear elastic 
stress analysis, no iteration will take place when this option is flagged. In 
inelastic analyses, the material tangent is interpolated and multiplied by the 
integration point strain which is directly sampled at the quadrature point. 
This option cannot be used for the advanced constitutive model since the cor- 
rect material tangent is not generated. 

*DISTRIBUTELOAD 

Body force and surface traction loadings are referred to as distributed loads 
in the MHOST program. The body force option includes gravity acceleration de- 
finable in any direction and centrifugal loading with the centerline and angu- 
lar velocity specified by the user. 

♦DUPLICATENODE 

The continuity of stresses at nodal points can be broken by defining two nodal 
points at the same geometrical location and connecting them by this option 
which enforces compatibility of displacements only. This option is used to de- 
fine the connections between generic modeling regions. 

♦DYNAMIC 

The generalized Newmark time integration algorithm is entered by setting this 
flag. A simple adaptive time stepping algorithm is employed at the user's op- 
tion. 
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♦ELEMENTS 

The elements Included in this version of the code are described in Table 
3.3-1 II. Core allocation is performed for the nodal and element quantities on 
the basis of maximum storage space requirements among the types of elements 
specified in this option. 

♦EMBED 

The subelement iteration capability is flagged by this option which signals 
the code to allocate the working storage for the subelement data in a hier- 
archical manner. The actual subelement mesh definition and the nodal and ele- 
ment data storage allocation take place when the individual subelements are 
defined. 

♦FORCES 

Concentrated nodal forces are defined and stored in an incremental manner. 
Core allocation takes place only when this option is invoked. 

*F RONT AL SOL UT I ON 

The frontal solution option for quasi-static analysis is implemented in this 
version of the code. Out-of-core storage devices are utilized and, hence, the 
capacity of the program is increased significantly. A direct access, rather 
than a sequential access, file is used for the solution buffer to avoid over- 
head due to the backspace operation in the back substitution phase. 

♦GMRS 

Generic modeling regions are defined as collections of elements that model 
geometrically parametrized parts of hot section components. Multiple generic 
modeling regions in a given mesh are connected using the duplicated node op- 
tion. Different parameters are specified for each generic modeling region, and 
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the Input data can be prepared separately. Internally, the complex of the ge- 
neric modeling regions is treated as a single mesh for the purposes of con- 
structing and solving the finite element equations. A table is prepared to 
report results separately for each generic modeling region. 

*L0U8IGNAC 

Parameters for numerical quadrature used in the mixed iterative processes are 
defined in a very precise way. Full integration, selective integration, or 
selective integration with filtering can be chosen for construction of the 
stiffness matrix. For residual vector integration, full and reduced integra- 
tion can be selected. The strain integration can be performed either by using 
uniformly reduced integration, trapezoidal integration with the reduced shear 
strain approximation or the previous quadrature with the filtering option. 

*M0DAL 

The free vibration modes of linear elastic structures are extracted when this 
option is invoked. The subspace iteration technique is utilized and a power 
shift option is included. 

♦NODES 

All the variables are defined and reported at nodal points. In the incremental 
processes, deformation and stress histories are integrated and stored only at 
the nodal points. Note that this architecture economizes storage substantially 
compared with fully integrated finite element displacement methods. 

♦OPTIMIZE 

A band width optimization procedure, based on the Cuthil -McGee algorithm, is 
applicable to in-core solution processes. No optimization procedure is re- 
quired for the frontal solver. 


C-5 L 
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♦PERIODICLOADING 


For transient calculations, nodal displacements and concentrated forces can be 
Input as sum* sol dal functions using this option. 

♦POST 

The postprocessing data, which contain all the information supplied to and 
generated by the code, are written to the file connected to FORTRAN unit num- 
ber 19. This file is formatted and can easily be manipulated by commercially 
available postprocessing packages with minor modifications. The header record 
of the file is designed to be compatible with the MARC post file which is pro- 
cessed by many finite element graphics packages. 

♦PRINTSETS 

The report generation is carried out on a nodal point basis with element inte- 
gration point options provided by interpolation using the shape function. 

♦REPORT 

The frequency of the line printer report generation is now controlled by this 
option. The default is to print at every increment. 

♦SCHEME 

Parameters that control the characteristics of the time integration operator 
are defined by the user. The default is the average acceleration algorithm 
commonly used in nonlinear dynamic finite element analysis. 

♦STRESS 

Boundary conditions for stress can be specified by the user as an option, al- 
though no mathematical justification is yet available for this type of con- 
straint. Any stress component can be prescribed at any nodal point. Simple 
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numerical tests have shown that inconsistent imposition of stress boundary 
conditions can lead to rapid divergence in the iterative process, so this op- 
tion should be used with caution. 

♦TANGENT 

This option invokes the modified Newton iteration procedure. The tangent ma- 
trix is updated until a user- specified iteration count greater than or equal 
to 1 is met. No updating occurs in subsequent iterations. The default itera- 
tion count is equal to 1, and the procedure generated by this value is known 
as the Kyi method of modified Newton iteration. This option has no effect 
when the BFGS process is employed. 

♦TEMPERATURE 

Nodal temperatures are read and used to generate thermal strains. These quan- 
tities are used also for the evaluation of creep strain and the Integration of 
coupled creep-plasticity models such as Walker's model. 

♦THERMAL 

Temperature dependent material properties are evaluated when this option is 
invoked and the appropriate user subroutine is provided to the system prior to 
execution. This operation is not necessary for the conventional creep model 
since temperature dependence can easily be incorporated into the user sub- 
routine CRPLAW. 

♦TRANSFORMATIONS 

Coordinate transformations at nodal points are specified by this option In 
which the angle of rotation is input by the user so that the code can generate 
necessary transformation matrices. The postprocessing file does not support 
coordinate transformations at this time. 
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♦TYING 


Multiple degree- of- freedom constraint equations are specified by the user 
through this option. Note that constraint equations are generated only for the 
displacement degrees-of- freedom. This option Is more flexible than the dupli- 
cate node option since constraints may be applied to individual degrees-of- 
freedom. 

The following parameters are used to signal to the MHOST program the presence 
of specific user subroutines: 

♦UBOUN *UF0RCE *UTEMP 

♦UCOEF *UH00K *UTHERM 

♦UUERIV *UPRESS 

As summarized above, the analysis capability included in the MHOST code covers 
most of the needs for inelastic analysis of turbine engine hot section compo- 
nents. The free format data input routines and report generation packages have 
been improved in the Task II program development effort to create a more com- 
fortable environment for users. Presently, the code consists of around 30,000 
lines of FORTRAN statements including extensive self-explanatory comment lines 
in each subroutine. The preliminary system document was generated in an auto- 
matic manner from these comment lines and the cross reference capability 
available on the Prime FORTRAN compiler and linker. 

Table 3. 3-1 1 1 summarizes the elements currently available in the MHOST code, 
Version 2.0, and the parameters used internally to define element character- 
istics. Further details are available in the MHOST User s Manual. 

Lagrangian rather than Serendipity quadratic functions are being employed in 
the MHOST code at present. This choice was based, in part, on the fact that 
Serendipity elements are known to behave in a less accurate manner than 
Lagrangian elements for nonlinear fracture mechanics applications based on the 
quarter point technique. In addition, the lumped mass approximation for nodal 
stress/strain recovery has not been implemented in Serendipity elements and, 
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hence, the accuracy and stability of these elements a 
context of mixed iterative computations. Lagrangian 


re not yet known in the 
elements have, in con- 


trast, already produced encouraging results in the mixed iterative environment. 


Table 3.3-III 

MHOST Code, Version 2.0 
Elements and Parameters 



Beam 

P. Stress 

P. Strain 

Axsym. 

Brick 

Shell 

ITYPE 

9 

3 / 101 

11/ 102 

10/ 103 

7 

75 

NELCRD 

3 

2 

2 

2 

3 

3 

NELNFR 

3 

2 

2 

2 

3 

6 

NELNOD 

2 

4/9 

4 / 9 

4/9 

8 

4 

NELSTR 

1 

3 

4 

4 

6 

8 

NELCHR 

3 

5 

5 

5 

5 

5 

NELINT 

2 

4 

4 

4 

8 

4 

NELLV 

3 

3 

3 

3 

3 

4 

NELLAY 

0 

1 

1 

1 

1 

5 

NDI 

1 

2 

3 

3 

3 

2 

NSHEAR 

0 

1 

1 

1 

3 

1 

JLAW 

1 

2 

3 

4 

5 

6 


NELCRD Number of coordinate data per node. 

NELNFR Number of degrees- of- freedom per node. 

NELNOD Number of nodes per element. 

NELSTR Number of stress and strain components per node. 

NELCHR Number of material property data for the element. 

NELINT Number of 'full' integration points per element. 

NELLV Number of distributed load types per element. 

NELLAY Number of layers of integration through the thickness of the shell 

element. 

NDI Number of direct stress components. 

NSHEAR Number of shear stress components. 

JLAW Type of the constitutive equation. 


Four-point Gaussian quadrature is employed in the stiffness calculations for 
the Lagrangian elements in accordance with common practice. The lumped mass 
matrix with unit density constructed for strain recovery is also based on re- 
duced four-point quadrature instead of the nodal quadrature used in the linear 
quadrilaterals. Similarly, other integrations in the strain recovery process 
employ reduced integration in order to avoid stress/strain oscillation. 
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The ongoing investigation of quadratic element technology will focus on repro- 
duction of nonlinear material behavior and sensitivity to Isoparametric dis- 
tortion. 

The algorithm used in the MHOST code to construct the residual load vector, 
which drives the iterative solution, has employed mixed strain/stress inter- 
polation for all the generalized strain/stress components. However, the mixed 
continuous strain recovery for certain terms can possibly violate the neces- 
sary condition for stability and cause convergence to an inferior numerical 
solution after a number of iterations. This numerical instability is particu- 
larly significant when the method is applied to incompressible problems and 
the dilitation terms are smoothed during the iterative displacement update. 
Due to the similarity in the mathematical structure of constraints, the trans- 
verse shear terms in the bending of plates and shells could cause the same 
difficulty if they are smoothed in the iterative process. Hence, a modified 
algorithm which utilizes the smoothed stress field for all the terms except 
the transverse shear so as to avoid possible numerical instability has been 
implemented into the MHOST code. 

3. 3. 5. 2 Data Storage and Control Structure 

The input data and information associated with the global finite element for- 
mulation and solution are stored in blank common using a conventional dynamic 
core allocation scheme. No array is prepared for quantities defined at element 
integration points. The storage required for the subelement iteration is allo- 
cated in the same work area in an indirect way. We prepare the element array 
for pointers and the actual work area is allocated when the storage is re- 
quired. At this time, the pointer array is filled. The conceptual representa- 
tion of the indirect core storage scheme enables us to set up a multi-level 
subelement refinement. 

In the recovery phase for residual vectors associated with the iterative solu- 
tion, a flag is set to indicate whether the global displacement mesh is sub- 
elemented or not. This flag is built into the code internally as part of the 
element definition data array. The same flag can also be attached to the de- 
finition data for subelements. From the global displacement mesh definition 
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upwards, the generic 
same fashion. That 
constructed. 


model i ng 
is, a 


regions (GMRs) can be handled internally in 
tree structure of element mesh topology 


the 

is 


The use of out-of-core storage may be needed to swap the information related 
to different mesh levels, and the storage scheme for the element level data 
manipulation is retained as it is. This means that the subelement mesh data 
are loaded to the array where the global mesh is stored when the subelement 
level analysis module is entered. If further refinement is activated, the 
software stack is pushed down again. Such a recursive storage scheme enables 
us to use the existing element manipulation scheme developed for the global 
solution without major modifications. 

The reporting scheme is such that the global solution is reported as it is 
carried out in the present code, and then the same data stack is used for 
generating reports separately at the subelement level. 

The nodal coordinates of the subelement mesh need to be stored in terms of the 
local coordinates defined in the global system. For catalogued mesh subdivi- 
sions such as the one shown on figure 3.3-4, these values are easily defined 
in the required system of coordinates. 

When the mesh is prepared by the user in physical space, the nonlinear equa- 
tion which represents its isoparametric counterpart is solved by an iterative 
manner. It is feasible to utilize the bilinear mapping function for the geo- 
metrical definition of subelements, assuming that in practice, the subelement 
size is small enough compared with the geometrical features of the entire 
structure. 


The actual subelement solution is carried out in the element residual calcula- 
tion routine. When the subelement flag is set in the element loop, the subele- 
ment solution driver routine is entered and the mixed iterative solution is 
performed in the subelement domain. A set of procedures, which are somewhat 
different from the global solution subroutines, has been developed to manipu- 
late matrices for the subelement domain. 
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3.3.6 Validation/Verification Analyses 


3.3.6. 1 Analysis of Plate and Shell Structures 

The complexity of plate and shell problems is mainly due to the higher order 
derivatives in the variational functional, or equivalently the necessity to 
introduce additional rotational degrees-of-freedom. The results of a series of 
analyses which demonstrate the characteristics of the MHOST shell element are 
presented in this section. 

Clamped Plate with Central Load 

The displacement convergence characteristics of the current version of the 
plate and shell element are illustrated on Figure 3.3-5. The solution for a 
clamped plate subjected to a concentrated load is plotted together with some 
well-known plate solutions obtained using the finite element method. In terms 
of the deflection at the center of the plate, mixed iteration produces results 
consistent with a conventional mixed formulation. Note that the Reissner- 
Mindlin plate solution (the result from displacement preconditioning ) provides 
a lower bound where the mixed and mixed iterative solutions provide an upper 
bound of the solution. The MHOST code could be tuned to reproduce the lower 
bound solution, if desired, by selectively undoing the continuous stress ap- 
proximation, particularly for the bending moments. 

Plate Supported at Corners 

The performance of the hourglass control algorithm is demonstrated using an 
example of a square plate supported at its corners and subjected to a concen- 
trated load at its center. This is the same model problem used by Belytschko 
and T say (Reference 12). The problem statement is illustrated on figure 3.3-6, 
and the results with and without the hourglass control are plotted on Figures 
3.3-7 and 3.3-8. The participation factor for the fully integrated transverse 
shear stiffness matrix was determined using this example and further valida- 
tion analyses were performed to ensure that the hourglass control modification 
does not significantly change numerical solutions. 
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Figure 3.3-5 Clamped Square Plate 
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Figure 3.3-6 Corner Supported Plate with Central Concentrated Load 



Figure 3.3-7 Lateral Displacement with Hourglass Control 
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figure 3.3-8 Lateral Displacement without Hourglass Control 
Simplified Turbine Vane Airfoil 


A highly simplified model of a turbine vane airfoil has been studied as part 
of the shell element validation effort. The finite element model and applied 
loads are shown on Figure 3.3-9. A primary motivation for analyzing this 
structure was to determine the effects of kink angles between adjacent ele- 
ments on the behavior of the MHOST mixed iterative solution process. As shown 
on Figure 3.3-9, kink angles as large as 40 “ were present in the model. 

Constant thickness analyses were performed using MARC element 75 (bilinear 
Reisner-Mindlin shell, Reference 11) and the NASTRAN QUAD4 element as well as 
the MHOST shell element. Typical displacement results at nodes 10 and 20 are 
shown in Table 3.3-IV. MHOST results with zero iterations compared very favor- 
ably with the MARC and NASTRAN values. The MHOST displacements for the iter- 
ated analysis were larger than the other tabulated results, but this behavior 
is not considered unreasonable for such a coarse mesh. However, the relatively 
slow rate of convergence of the iterative process (14 iterations) is a source 
of concern. 
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Table 3.3-IV 

Displacements for Constant Thickness Vane Model 


Deck/Element 

Number of 
Iterations 

x-Displacement (in.) 

Node 10 Node 20 

MARC/75 

- 

0.469 x 10- 2 

0.745 x 10-2 

NASTRAN/QUAD4 

- 

0.477 x 10-2 

0.764 x 10-2 

MHOST/75 

0 

0.468 x 10-2 

0.762 x 10-2 

MHOST/ 75 

14 

0.496 x IQ" 2 

0.822 x 10-2 


The investigation was continued with analyses of the variable thickness model. 
In this case, results using MARC element type 4 (Hermitian curved shell ele- 
ment) were obtained in place of QUAD4 results. Displacement solutions obtained 
for nodes 10 and 20 are shown in Table 3.3-V. The tabulated variable thickness 
results show the same trends as those exhibited by the constant thickness dis- 
placements. In an effort to improve the convergence results, continuity of 
stresses/strains across the largest kink angle was destroyed by introducing 
duplicate nodes coincident with nodes 6-10. This action had the effect of re- 
ducing residuals by about an order of magnitude as is shown on Figure 3.3-10. 

_2 

The displacements in the x direction at nodes 10 and 20, 0.320 x 10 inch 
and 0.400 x 10”^ inch respectively, compared favorably with corresponding 
results for the continuous stress/strain case (Table 3.3-V). 

In summary, the cantilevered vane results showed that, while substantial kink 
angles between adjacent shell elements do not cause the MHOST iterative pro- 
cedure to produce erroneous results, such angles can slow convergence to the 
point that complete stress/strain continuity constraints may have to be re- 
laxed on a selective basis in order to achieve acceptable convergence rates. 
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Table 3.3-V 

Displacements for Variable Thickness Vane Model 


Deck/El ement 

MARC/4 

MARC/75 

MHOST/75 

MHOST/75 


Number of 

x-Dlsplacement (Inch) 

Iterations 

Node 10 

Node 20 

- 

0.267 x 10- 2 

0.356 x 10-2 

- 

0.274 x 10“2 

0.357 x 10-2 

0 

0.266 x 10-2 

0.341 x 10-2 

20 

0.331 x 10-2 

0.415 x 10-2 



Figure 3.3-10 Convergence of Residual Forces 
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Transient Response of a Simply Supported Bar 


The pinned beam shown on Figure 3.3-11 represents a simple validation case for 
the MHOST time integration algorithm (constant average acceleration option). 
The beam was modeled using three MHOST shell elements, and was subjected to 
the impulsive pressure load shown on the figure. The lateral displacement time 
history is compared to a MARC modal solution (Timoshenko beam element model) 
on Figure 3.3-12. As the plotted results show, the MHOST solution is in excel- 
lent agreement with the MARC results over the complete time range (0.0 to 0.18 
seconds). 
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Figure 3.3-11 Model and Dynamic Loads for Pinned Beam 
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Free Vibrations of a Cylindrical Fan Blade 


The MHOST free vibrations capability was verified by analyzing the well-known 
Lindberg and Olson cylindrical fan blade (Reference 49). The MHOST shell ele- 
ment model of the fan blade as well as a NASTRAN model that employs an 8 node 
element (Reference 50) are defined on Figure 3.3-13. The MHOST frequencies for 
the first ten modes are compared to NASTRAN and test values in Table 3.3-VI. 
Overall, the MHOST frequencies are in good agreement with test values, and re- 
present approximately the same level of approximation as the tabulated NASTRAN 
results. Macroscopic comparisons between MHOST and experimental mode shapes 
are shown on Figure 3.3-14. The MHOST program captures the essential charac- 
teristics of all the modes shown on the figure. 
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figure 3.3-13 Models of Lindberg and Olson Fan Blade 
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Table 3.3-VI 

Frequencies for Lindberg and Olson Fan Blade 




Frequencies (Hz) 

lode 

iml 

Test 

NASTRAN 

1 

A 

85.6 

86.7 

2 

s 

134.5 

139.6 

3 

S 

259.0 

250.3 

4 

A 

351.0 

349.6 

5 

S 

395.0 

407.9 

6 

A 

531.0 

552.5 

7 

A 

743.0 

777.2 

8 

S 

751.0 

759.7 

9 

S 

790.0 

826.2 

10 

A 

809.0 

920.9 

*S = 

symmetric 



A = 

anti -symmetric 




V\ 


Figure 3.3-14 Comparisons Between MHOST and Experimental Mode Shapes 





3. 3. 6. 2 Analyses with Subelement Iterations 


Cantilever Beam 


The cantilever beam shown on Figure 3.3-15 was analyzed using a 2 x 10 global 
mesh in conjunction with uniform 2x2 subelement meshes. The subelement solu- 
tion extracted exact bending stress/strain and displacement values at all 
nodes including the interior subelement mesh points. 

A global 4 x 20 model of the beam was also analyzed to provide reference 
timing information. The computational efforts required for the standard and 
subelement approaches are summarized in Table 3.3-VII. It should be noted that 
the subelement procedure offers a computational advantage over the global ap- 
proach in large-scale applications even though the total computational efforts 
listed in Table 3.3-VII do not reflect this fact. Since matrix solution costs 
increase by powers of the number of degrees-of- freedom, the reduction in this 
cost obtained using the subelement approach in large-scale applications will 
more than compensate for the cost increase associated with subelement calcula- 
tions which varies linearly with problem size. The differences in matrix solu- 
tion costs associated with the subelement and global procedures are apparent 

even in the cantilever example as is shown in Table 3.3-VII. 



1000 LB 
1000 LB 


FIGURE NOT TO SCALE 


E _ 10 6 psl • SUBELEMENT MESH 

u-0.3 
t - 1.0 INCH 

CANTILEVER BEAM PROBLEM (PLANE STRESS) 

Figure 3.3-15 Cantilever Beam Validation Problem 
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Table 3.3-VII 


Central Processing Unit (CPU) Time for the 
Beam Validation Problem 


Run A 2x10 mesh with 2x2 uniform subelement division throughout 66 main 
degrees-of- freedom and 360 subelement degrees-of- freedom 

Run B 4x20 mesh of the standard finite elements 210 degrees- of- freedom 


CPU Time (seconds In PRIME 9950) 


Operation Run A Run B 

Matrix Assembly 0.803 2.267 

Matrix Solution 0.091 2.290 

Residual Calculation 10 . 803* 3.906 

TOTAL 11.607 8.463 


♦Including the subelement solution which averages 0.5 sec/global 
element. 

Plane Stress Problem with a Reentrant Corner 

The plane stress problem with a reentrant corner shown on Figure 3.3-16 was 
analyzed as a validation case for problems with discontinuities. Related work 
at Swansea has indicated that stress fields obtained from the mixed iterative 
procedure are oscillatory when stress singularities are present. These numer- 
ical difficulties are related to the stability condition of Babuska-Brezzi , 
and the possible unstable characteristics of equal order interpolation for 
displacement and stress have been pointed out by Oden (Reference 33). The 
stress distributions obtained by MHOST after global iteration with continuous 
stress fields are shown on Figures 3.3-17 through 3.3-19. Except for a small 
amplitude stress oscillation, no significant indication of numerical insta- 
bility is observed. 


The nature of the discontinuity at the reentrant corner is such that the con- 
tinuous stress constraint is excessive. When this constraint is removed, a 
more accurate representation of singular behavior at the corner is obtained as 
shown on Figure 3.3-20. A further refinement shown on Figure 3.3-21 is ob- 
tained by using a 2 x 2 subelement mesh representation in the three elements 

at the singular point. 
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figure 3.3-17 Stress Distribution (o x ) for L-Shape Domain with Continuous 
Stress field 
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Inelastic Analysis of a Plate with a Hole 


An elastic-plastic analysis of a plate with a circular hole was performed 
using the model shown on Figure 3.3-22. Elements 7, 8, 17, and 18 were divided 
into 2x2 uniform subelement meshes. The material data used in the analysis 
are shown in Table 3 . 3— VI II. In this problem, an average of 4 to 5 iterations 
was required in each subelement to reach convergence with a displacement 
tolerance of 0.1. The deformed shape after one increment is plotted on Figure 
3.3-23. As shown in Table 3.3-IX, the equivalent plastic strain detected at 
the subelement level in element 8 was an order of magnitude larger than the 
value observed in the global solution. This result indicates the potential 
capability of the subelement approach in nonlinear fracture mechanics appli- 
cations. 



Figure 3.3-22 Finite Element Model of a Plate with a Circular Hole 
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Table 3.3-VIII 


Parameters for the Validation Model 
of a Hole In Finite Plate 


Thickness 
Young's Modulus 
Poisson's Ratio 
Yield Function 
Equivalent Stress (psi) 

3.0 x 10 4 
3.5 x 10 4 

4.0 x 10 4 

5.0 x 10 4 

6.0 x 10 4 


1.0 inch 

3.0 x 10' psi 
0.3 

Plastic 

Equivalent Strain 


u • V 

0.5 x 10-2 
1.5 x 10-2 
1.0 x 10- 1 
1.0 



Figure 3.3-23 Deformation After One Increment 
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Table 3.3-IX 


Solution at Point A by the Global and 
Subelement Iterations After Increment One 
Total Load F * 1.1 F 0 


x Displacement (inches) 
Von Mises Stress (psi) 

<* x (psi) 

(psi) 

T xy <psf) 

Equivalent Plastic Strain 


Global Node 
6.5360 x 10~ 4 
3.0003 x 10 4 
-1.9666 x 10 3 
-3.0114 x 10 4 
4.0285 x 10 3 
3.5822 x 10' 6 


Local Node 
6.5360 x 10“ 4 
3.0185 x 10 4 
-2.7765 x 10 3 
-3.0941 x 10 4 
3.2650 x 10 3 
1.9407 x 10‘ 4 


3.3.7 Three-Dimensional Sol id Element Analyses 
Burner Blister Specimen 

A series of elastic-plastic analyses has been completed for a burner blister 

specimen configuration. The blister specimen was modeled as a 45° in-plane 

wedge with a radius of 1.5 inches. Figures 3.3-24 and 3.3-25. The thickness 
for the out-of-plane direction was 0.05 inch and three planes of nodes (each 
plane containing 199 nodes) were used through the thickness. The three planes 
were designated as bottom plane, mid-plane, and top plane respectively. Roller 
in-plane boundary conditions were applied along the horizontal and 45“ wedge 
boundaries. Out-of-plane boundary conditions prevented rigid body motion. 
Eight-noded three-dimensional brick elements were selected for both the F410ST 
and MARC analyses, with two elements used through the out-of-plane thickness 
direction. Elastic and plastic material properties corresponded to realistic 
burner liner materials. The external loading involved a radially varying tem- 
perature distribution, which ranged from 1800“F at the vertex hot spot to 

1100“F at the outer radius, Figure 3.3-26. The distribution was based upon 
experimental burner temperature measurements, was applied proportionally for 
incremental loading, and was linear through the thickness direction. 
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figure 3.3-24 Burner Blister Specimen Model In-Plane General Breakup (Three 
Dimensional ) 


VERTEX 


0.025 


0.1 



figure 3.3-25 Burner Blister Specimen Model In-Plane Local Vertex-Neighbor- 
hood Breakup (Three-Dimensional) 




RADIAL LOCATION 

figure 3.3-26 Hot Spot Temperature Distribution (Mid-Plane) for Burner Blis- 
ter Specimen Model 

In addition to a three-dimensional modeling of the burner blister specimen, a 
finer-meshed axisymmetric blister specimen model. Figure 3.3-27, was also de- 
veloped and run on MARC. The axisymmetric model contained four four-noded ele- 
ments through the thickness direction, as well as a mesh breakup near the ver- 
tex that was somewhat finer than the 3-D model mesh breakup near the vertex 
(compare Figures 3.3-27 and 3.3-25). The results associated with the MARC axi- 
symmetric model calculations were taken to be base case numbers, to be used to 
judge the accuracy of the MHOST and MARC three-dimensional model calculations. 
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Figure 3.3-27 Axi symmetric Burner Blister Model for MARC-Breakup Near Vertex 

Agreement of results between the MHOST and MARC three-dimensional model calcu- 
lations was very good for most quantities. Typical of this good agreement are 
the results shown on Figure 3.3-28, where the Mises plastic strain at the mid- 
plane vertex node is plotted from both sets of results. As can be seen from 
this figure, the comparison is excellent even at 100 percent load for which 
significant plastic strain occurs at the vertex node. Even better agreement 
was obtained for the in-plane radial displacements which sometimes showed 
three or four figure exact matchup between MHOST and MARC values. It should 
also be noted that whenever good agreement occurred between MHOST and MARC 
three-dimensional calculations, the corresponding quantity found from the MARC 
axi symmetric calculation also agreed very well with both 3-D results. 

On the other hand, poor agreement between the MARC and MHOST three-dimensional 
calculations was found for the out-of-plane displacement values, in particular 
for the blister displacement at the vertex nodes as indicated in Table 3.3-X. 
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Figure 3.3-28 Mises Plastic Strain at Mid-Plane Vertex Node 

Table 3.3-X 

Out-of-Plane Displacement at Mid-Plane Vertex 
Node for 100 Percent Loading 

MARC AXISYMMETRIC: 5.0982 x 10“ 3 

MARC 3-D: 1.2473 x 10“ 3 

MHOST 3-D: 4.7697 x 10‘ 3 


In general, the MHOST 3-D out-of-plane displacement values were about three or 
four times greater than those given by the MARC 3-D calculations, but were In 
quite good agreement with the MARC axl symmetric calculations, as Indicated In 
Table 3.3-X. In fact, results from the MARC axlsymmetrlc case showed good to 
excellent agreement with all correspondl ng MOST 3-D results for all quanti- 
ties. 

It was concluded, therefore, that the MHOST code gave good to excellent re- 
sults for all quantities in the three-dimensional burner blister specimen 
model with two three-dimensional elements through the thickness. In contrast, 
two three-dimensional elements through the thickness was not a fine enough 
breakup for MARC to generate accurate out-of-plane displacement values. It is 
felt that a possible cause for this poor accuracy is that the element coarse- 
ness through the thickness resulted in too high an out-of-plane stiffness for 
the MARC calculational procedure to be able to handle properly. It appears, 
therefore, that the MHOST code can handle coarser meshes better than the MARC 
code. The primary reasons for this advantage can be found from among the fol- 
lowing factors: integrated use of reduced integration element technology, 
nodal evaluation of constitutive relationships, and equilibrium (Loubignac) 
Iteration. 

With regard to computer running times required for MHOST and MARC 3-D model 
calculations, MHOST appears to be somewhat faster than MARC for identical 
problems, geometries, loadings, etc. Time comparisons between the two programs 
are shown in figure 3.3-29, where 100 percent of the load was applied in one 
step (increment) and the number of iterations was varied. The plots here show 
the MHOST code to be about 40 percent faster than the MARC code for this par- 
ticular problem setup and loading. This type of time savings for MHOST over 
MARC (i.e., about 40 percent) is typical of the comparisons that were made, 
regardless of whether the loading involved one increment and a number of iter- 
ations, or vice-versa, a number of increments (adding up to 100 percent load) 
and zero or one iteration per increment. 
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Figure 3.3-29 Computer Running Times for 100 Percent Load in One Step (Three- 
Dimensional Model) 

from the discussion presented in the preceding paragraphs, it can be concluded 
that the MHOST code, when compared against the MARC code, provides significant 
accuracy improvement at reduced cost for applications of the blister specimen 
type. 

NASA Benchmark Notch Specimen 

Three-dimensional elastic-plastic analyses of the NASA Benchmark Notch Speci- 
men have been performed using the MHOST code. In order to avoid duplication, 
the finite element results for this specimen are discussed together with 
boundary element method results in Section 3.4. 5. 7. 
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3.3.8 List of Symbols 


Alphabetical 

Symbol 

a 

*"w 

a j 

B 

rw 

c 

C 

0 

°ijki 

e 

E 

£ R 

X 

i 

fj 

G 

#"W 

g 

**w 

h 

I 

J 

K 

rw 

L 

M 

rw 

N 

rw 


Description 
Acceleration vector 
Acceleration vector component 
Strain-displacement matrix 
Damping matrix 
Geometric constant 
Material modulus matrix 

Fourth order tensor component of material modulus 
Nodal strain vector 
Young's modulus 

Prescribed surface traction vector 

Body force vector 

Node point force vector 

Body force vector component 

Gramm matrix for the finite element basis 

Gravity acceleration vector 

Thickness for plates and shells 

Functi onal 

Jacobian matrix for the isoparametric transformation 
Displacement stiffness matrix 
Lagranglan functional 
Mass matrix 

Finite element basis function 
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List of Symbols (continued) 


Alphabetical 

Symbol 

n 

"j 

P 

Q 

Q' 

R 

R 

$ 

/■w 

r 

t 

u 

U 

u j 

V 

v 

ja 

v j 

w 

w 

X 

x j 


Description 
Unit normal vector 

Cartesian component of the unit normal 
Penalty matrix 
Initial strain terms 

Nodal forces generated by initial strains 
Radius 

Rotation tensor 
Nodal stress vector 

Radial coordinate for the axi symmetric geometry 

Traction vector 

Traction vector component 

Right stretch tensor 

Displacement vector 

Displacement vector component 

Space for admissible displacement variation 

Velocity vector 

Vector for BFGS update 

Velocity vector component 

Lateral deflection 

Vector for BFGS update 

Position vector 

Cartesian component of the position vector 


List of Symbols (continued) 


Greek Symbols 

a,/3,T 

l 

6 

e 1j 

£ * n * i 

e 

X, u 

x 


V 


p 


°1j 

t. • 
fc lj 

d 

S2 


3J2 


Description 

Constants 

Kronecker's delta 

Displacement at a given point 

Infinitesimal strain tensor 

Isoparametric coordinates 

Constant 

Lame-Navi er constants for the isotropic elasticity 

Load factor for the arc-length method 

Poisson's ratio 

Material density 

Stress tensor 

Stress tensor component 

Deviatoric stress tensor 

Deviatoric stress tensor component 

Null set 

Problem domain 

Boundary of the domai n 


Sub- and 
Superscripts 

e 

h 

i » j *k » • • • 


I ,J ,K, . . . 


Description 

Quantities defined on each element 
Finite element approximation 

Indices for vector and tensor component (when used as 
subscripts); iteration and incrementation counter (when 
used as superscripts) 

Nodal point counter 
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List of Symbols (continued) 


Sub- and 
Superscripts 

L 

MAX 

s 

S 

T (Sub) 

T (Super) 


Description 

Local subelement mesh 

Maximum number 

Shear and transverse shear component 

Subel ement 

Tangent 

Transpose 

Initial quantities 


0 
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3.4 BOUNDARY ELEMENT METHOD 


3.4.1 Overview 

It Is the goal of the advanced formulation development portion of this program 
to develop a computational technique for the solution of linear, nonlinear and 
transient problems in gas turbine engine hot section components. This tech- 
nique is to be distinct from, and complementary to, the Finite Element Method. 
The existence of such a computational tool will enhance the ability to cali- 
brate the other codes developed under this contract. In addition, it is to be 
expected that different techniques will prove optimal, in terms of efficiency 
or accuracy, for particular types of component analysis. Since almost all gen- 
eral purpose structural analysis computer programs presently available employ 
the displacement finite element method, the new program developed as part of 
the advanced formulation development effort can be expected to extend the 
ability to perform realistic analyses of hot section components. 

During the first year of this program (Task IC), Pratt & Whitney and its sub- 
contractor, the State University of New York at Buffalo (SUNY-B), developed a 
new general purpose structural analysis program, BEST3D (Boundary Element 
Stress Technology), based on the use of the Boundary Element Method (BEM). 
During this work, the boundary element method was implemented for very general 
three-dimensional geometries, for elastic, inelastic and dynamic stress analy- 
sis problems. 

In the second year of the program (Task IIC), Pratt & Whitney and SUNY-B have 
continued the theoretical and numerical development and the computer implemen- 
tation of the BEM, making very significant advances in a variety of areas. 
Major developments accomplished during Task IIC include: 
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1. Incorporation of substructuring capability in the nonlinear BEM 
stress analysis. 

2. Algorithm development, coding and validation of an embedded time 
algorithm for elastodynamic problems. 

3. Development of a method for representing the effects of embedded dis- 
continuities without explicit boundary modeling. 

4. Development and two-dimensional testing of an algorithm for the cal- 
culation of the natural frequencies and mode shapes of a structure, 
and the preliminary design for incorporating this capability in 
BEST3D. 

5. Significant improvement in the accuracy and efficiency of many of the 
numerical algorithms in BEST3D. 

6. Verification of the nonlinear solution capabilities of BEST3D using 
externally generated data. 

The second year development effort is discussed in more detail below. Section 
3.4.2 updates the results of the BEM literature survey originally conducted 
during Task IC. Section 3.4.3 summarizes new developments in the analytical 
and numerical formulation of the BEM for elastic, inelastic and dynamic prob- 
lems in three dimensions. Modifications to the BEST3D code are discussed in 
Section 3.4.4. Validation/verification of BEST3D is discussed in Section 3.4.5. 
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3.4.2 Literature Survey Update 


During 1984, no boundary element papers have been published that are likely to 
cause significant redirection of work In the Inelastic Methods Program. It 
appears that the boundary element technology used In BEST3D for multiregion 
Isotropic and anisotropic elastic stress analyses, multireglon dynamic 
(steady-state and transient time domain) analyses and multiregion Inelastic 
analyses under monotonic and cyclic loading remains at least current with the 
published literature. 

A comprehensive textbook by Brebbia et al (Reference 1) and an advanced mono- 
graph edited by Banerjee and Mukherjee (Reference 2) appeared in 1984, and 
essentially document the state-of-the-art. Additionally, Ingrafia (Reference 
3) described the use of special shape functions for fracture mechanics analy- 
sis, along lines originally explored by Cruse and Wilson (Reference 4). Rizzo 
et al (Reference 5) described dynamic analysis of some earth structures prob- 
lems involving single homogeneous regions using the Fast Fourier Transform. 
Brebbia and Nardini (Reference 6) have explored the calculation of natural 
frequencies in single region, two-dimensional boundary element analysis. 

3.4.3 Formulation Development 

3. 4. 3.1 Summary 

Important advances have been made during Task IIC in extending the BEM formu- 
lation for three-dimensional stress analysis of gas turbine engine structures. 
The most significant formulation advances have been in the area of dynamic 
analysis, where a real variable, time embedded technique has been developed, 
in the calculation of natural frequencies and mode shapes, in the representa- 
tion of multiple embedded discontinuities and In the basic understanding of 
numerical algorithms employed in the BEM. These extended or newly developed 
formulations are discussed in the subsections below. The basic BEM formulation 
Is only very briefly reviewed, as full details are available in the First 
Annual Status Report (NASA CR-1 74700). 
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3. 4. 3. 2 Linear and Nonlinear Stress Analysis 


3. 4. 3. 2.1 Review 


By making use of the reciprocal work theorem, the governing differential equa- 
tions for a three-dimensional (homogeneous) structure under combined thermal, 
mechanical and body force loadings can be converted to an integral equation 
written on the surface of the structure. This integral equation is: 

c u u i - s tG u *i - F u u i ] ds * / f i **>y> d * f 3 - 4 - 1 * 

where Wj = 6^, y3 = coefficient of thermal expansion. 


T ikj - the stress »,-k 


due to a point force e., 

J 



Tikj and ^ij are defined in the F^st Annual Status Report 
(NASA CR- 174700), 


and 


ij 


5 *- D ijk t k " S ijk u k ] ds * / (T ijk f k + M ij T) dv 


(3.4-2) 


allows calculation of stresses at any interior point where they are required. 

A similar equation for interior displacements can be obtained by setting c^. 

= 6 .. in equation (3.4-1). 

^ J 

In a purely elastic problem BEM stress analysis can be carried out entirely on 
the boundary of the structure. Once a physically reasonable set of boundary 
conditions has been prescribed, equation (3.4-1) can, in principle, be solved 
for all of the remaining boundary displacements and tractions. 


It is generally impossible to solve equation (3.4-1) exactly for real struc- 
tures and loading conditions. Suitable approximations of the boundary geome- 
try, displacements and tractions must be used in order to reduce equation 
(3.4-1) to a system of algebraic equations. The present version of BEST3D 
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models boundary geometry and boundary values of field quantities using linear 
and/or quadratic Isoparametric shape functions. The surface Integrals In equa- 
tion (3.4-1) are then evaluated numerically using product Gaussian quadrature 
rules. The numerical implementation of the BEM Is discussed In detail In text- 
books by Banerjee and Butterfield, and by Brebbla as well as In the First 
Annual Status Report (NASA CR-174700). 


In the case of inelastic analysis, the volume Integrals In equation (3.4-1) 
cannot be calculated a priori , since they require knowledge of Inelastic 
strain, which is itself a part of the solution. In this case equations 
(3.4-1), (3.4-2) and the inelastic material model can be regarded as a coupled 
system of nonlinear equations. In the numerical implementation of the BEM 
equation (3.4-2) is used to calculate the stresses at interior points, and the 
nonlinear material model is then used to evaluate inelastic strain. Since the 
volume integrals of inelastic strain vanish except in regions of nonlinear 
material response, approximations of geometry and field quantities are re- 
quired only where nonlinearity is expected. In BEST3D volume models utilize 
quadratic isoparametric shape functions or simple rectilinear cells. In the 
first case the isoparametric shape functions are also used to model inelastic 
strain variation. In the second case an exact solution for the uniform initial 
strain problem, discussed in Section 3. 4. 3. 2. 6, is used. 

The remainder of Section 3. 4. 3. 2 discusses significant new developments 
carried out during Task IIC. 


3. 4. 3. 2. 2 Redefinition of Plasticity Formulation 


In the First Annual Status Report (NASA CR-174700) the plasticity formulation 
was defined with respect to a known distribution of initial strain. It has 
been found to be more convenient to define the displacement field in terms of 
the initial stresses defined over the volume, so that: 


'ij 


u 1 (x o> 


Jr 


S CG 1J (x ’ x o ) V x) ‘ 


F ij .(x,x Q )u 1 (x)dr 


V B 1kj (z *^ 5 ik (z)dv 


(3.4-3) 


where 


B ikj 


1 

1 6tt n 


n 




[(l-2v) U 1k — - « jk 7 


y 1 , , 3 Wk , 

-=-) - 6 ,• -i ’s — J 


ij r 
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The stresses at an Interior point can be obtained from 


’jk 


. / 

~ S 


^ D ijk t i 


f 

~ u.] ds + ;, E 


ijk i 


.0 


V ipjk ip 


dv 


+ 1 • o 

J ipjk “Ip 


(3.4-4) 


where the last term on the right-hand side of equation (3.4-4) is the princi- 
pal value. The kernels D and S are defined in Banerjee and Butterfield and 


E ipjk = ' ( ” 1)r « [c 4 (s 1p *jk - 5 ij s kp ■ *ik 4 jp ' “ 4 Jk Vp 1 


“ 4 1p */k 


av(6 ij Vp * 4 ik y j y p ’ { jp Vk * *kp y t yj> * VjV P ] 


C 3 = -1/4* (1-v) 

[(a^ - 2) - v(a^ - 4)] 6^ + [1 - v(o + 2) d^] 

= -1 / (1-v) , a = 3 

Equations (3.4-3) and (3.4-4) provide the formal basis for developing the 

plasticity algorithm. The initial stresses a?, defined in equations 

* 

(3.4-3) and (3.4-4) include the effects of all inelastic strains (thermal, 
plastic and creep). The magnitude of the initial strain, however, is not known 
a priori (except for the thermo-elastic problem) and must be determined by 
satisfying the constitutive relations. 

3. 4. 3. 2. 3 Plasticity Solution Sequence 

The displacement and stress equations at the boundary (equations 3.4-3 and 
3.4-4) and interior points can be assembled to arrive at the system equations 
as: 
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u b i b . T b u b ♦ E b i-Xh* 1 ' = 0 


= U 1 t b - T 1 u b * E 1 6 o(b,1) (3.4-5) 
4 = U a t b - T° u b + E° S 0<b * i) 

where u , t , a 0 are the displacement, traction and initial stress increment 
vectors. The superscripts b, i in the above equations indicate reference to 
boundary and interior node related equations. 

By collecting the known (y) and unknown (x) values of traction and displace- 
ment rates and their coefficients together, the above equations can be recast 
as: 


A b i = BS * C b o° = b b ♦ b 0b 

u 1 ' = A 1 x + B 1 y + C 1 o° = A 1 x + b 1 + b 01 (3.4-6) 


The algorithms described yield the solution of the system defined by (3.4-6) 
together with the constitutive model. This includes complete knowledge of the 
initial stress distribution o° within the yielded region. Since a 0 is not 
known a priori for a particular load increment and since the complete system 
is nonlinear, an iterative process is employed within each loading stage. 


An important feature incorporated in all the iterative algorithms in the pre- 
sent work is the utilization of the initial stresses generated by the past 
history. In this procedure, the path followed by the previous load increment 
is used to extrapolate the initial stresses at the beginning of the current 
increment before the iterative operations. This substantially reduces the com- 
puter time. 
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The incremental algorithm can be described as follows: 


(a) Obtain the elastic solution for an arbitrary increment of boundary 
loading y from 


i = [A b ] -1 b b 

u 1 = A 1 x + b 1 
and 

o = A a x + b° 

(b) Scale the elastic solution such that the highest stressed node is at 
yield. 

(c) Impose a small load increment y (usually about five percent of the 
yield load) and an initially estimated value b° based on the pre- 
viously generated history of initial stress and evaluate 


i . [A b ] -1 b b * b° b 
and 

b = A° x + b° + b 0a 

It should be noted that at the first load step, the initial stresses are 
zero since no prior plastic history exists. 

(d) Accumulate all incremental quantities computed in the previous step. 

(e) Evaluate the current constitutive matrix using the new history and 
incremental values and calculate the initial stress rates and accumu- 
late them. 
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(f) Return to step (c) and apply the next load step with non-zero b°» 
If the current Increments of Initial stresses are less than a suit- 
able error norm (normally 0.005 times the yield stress). 

(g) Return to step (c) and compute the Incremental quantities due to ini- 
tial stress only (i.e., boundary loading change is zero). 

If the number of iterations executed is greater than a specified limiting 
value (usually 50), the system is assumed to have reached the state of col- 
lapse. 

The plastic solution algorithm is summarized in Figure 3.4-1. 



Figure 3.4-1 Plasticity Solution Sequence 
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3, 4, 3, 2, 4 Inelastic Material Models 


The three constitutive models which are presently included in BEST3D are: 

(a) Von Mises model with isotropic variable hardening 

(b) Two- surface model for cyclic plasticity 

(c) Combined plasticity and creep model of Walker. 

Of these, the models (a) and (c) are described in the earlier sections of this 
report. The two-surface model used here is similar to that described by Krieg 
except that a new form of hardening function has been introduced. 

In the two-surface model, it is assumed that during loading in the plastic 
state the stress state remains in contact with the inner yield surface, known 
as the loading surface. The outer surface, known as the limiting or bounding 
surface, is allowed to follow hardening rules identical to those of elastic, 
linearly hardening isotropic or kinematic plastic theory. 

The hardening of the loading surface depends on a distance vector that joins 
the stress state to the bounding surface such that the system is very stiff 
when the stress state is remote from the bounding surface and assumes essen- 
tially the hardening values of the bounding surface when the stress state is 
in the proximity of this surface. This variation of stiffness from large to 
small asymptotic values allows a smooth transition while the stress state 
moves through three distinct zones: an elastic region which is associated with 
recoverable strains when the stress state is within the loading surface, a 
plastic region where the system adopts the stiffness of the bounding surface, 
and an in-between metaelastic region when the stress state is on the loading 
surface but inside the bounding surface. 
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The behavior In the elastic and plastic region Is governed by the stress 
strain relations: 

? e <&> <£*> T ? e 

do = D e FT'+TI' 

p e 

where o* is the reference point on the boundary surface corresponding to the 
stress point a, (as defined by Mroz in 1967). The function F used in this de- 
velopment is the standard yield function of Von Mises. 


de 


(3.4-7) 


In equation (3.4-7) the terms H p and H g are defined as: 


H e • <£*» 


<*F \T p*e j3F j 


ao’ 


and 


H = h(— ) n , h, n and o ref are material parameters 
P °ref 

a is the distance (in six-dimensional stress space) between the current 
stress point and the last elastic location of the loading surface 
center. 

Any unloading during the deformation is indicated by 
dl * n T do < 0 where n is the unit normal at a* 

n = (|Ey) / [(3F/3o*) T (3F/3 o*)] 1/2 


which leads to a purely elastic response until the loading path touches the 
loading surface when once again a new reference point a,* is established. 
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3. 4. 3. 2. 5 Multi-GMR (Generic Modeling Region) Plasticity Algorithm 

In the case of a multi region problem, the system matrices for the boundary 
nodes for each region are organized as: 

Ax + By + Co 0 = 0 


or Ax = -(By ♦ Co 0 ) = b (3.4-8) 

MM MM M 

Equation (3.4-8) is identical in nature to the corresponding system equations 
for the multi -GMR elastic problem and therefore can be used in the multi region 
assembly and solution scheme without any change. The iterative process used 
does require an efficient process for the repeated solution of the algebraic 
system (3.4-8). 

3. 4. 3. 2. 6 Embedded Hot Spots and Discontinuities 

An exact solution for the case of uniform initial stress or initial strain in 
a rectilinear cell of dimensions a x b x c has been derived, based on the 
stress solution in Reference 7. This solution can be expressed as: 

Uj (0 = B? kj (x,i) a° k (x) (3.4-9) 

and a j -| ( £ ) = M ijkl (x,?) o? k (x) (3.4-10) 

It is important to note the following features of the above equations for dis- 
placements and stresses: 

(a) there is no integration required 

(b) equations (3.4-9) and (3.4-10) can be added directly to the respec- 
tive boundary integral equations for the displacements at the bound- 
ary and interior points and for the stresses at interior points. 
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For a hot spot equations (3.4-9) and (3.4-10) simplify to: 

UjU) = H.(x,l) T(x) (3.4-11) 

a.jU) = J.j(x,|) T(x) (3.4-12) 

where T(x) Is the temperature of a hot spot of volume a x b x c. In the pre- 
sence of hot spots equation (3.4-8) Is modified to 

Ax + By + Co 0 + HT = 0 (3.4-13a) 


or Ax = -(By + Ca° + HT) = b 


(3.4-13b) 


Equation (3.4-13) can be formed for each generic modeling region (GMR) and 
used in the solution. Discontinuities are allowed without having any formal 
discretization by simply assuming that the stresses within the volume occupied 
by it are zero. Thus, for example, in an elastic system if we introduce an 
arbitrary initial stress system £° within the discontinuities, then the 
boundary displacement equations and the interior stress equations at a point 
within a discontinuity can be written as: 

Ax + By + Ca° = 0 (3.4-14a) 


a = 0 = A°X + B a y + Ma° 


(3.4-14b) 


By eliminating a 0 between equations (3.4-14a) and ( 3 .4- 14b) , 
Ax + By - CM" 1 (A°x + B a y) = 0 

or (A - CM'V) x = (CM -1 B a - B) y 


or Ax = b 


(3.4-15) 
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Equation (3.4-15) can be formed for each GMR and solved as before. 

It should be noted that displacements and stresses near the void change sig- 
nificantly as a result of its presence. This is due to the fact that both the 
boundary solution U,y) and the initial stress (a Q ) occur in both equations 
(3.4-14a) and (3.4-14b). 

3. 4. 3. 2. 7 Anisotropic Elastic Materials 

The ability to perform stress analysis for materials whose elastic behavior is 
anisotropic is of considerable importance in the study of hot section compo- 
nents. Such anisotropy can occur in three ways: 

(1) Because processing (such as rolling or forging) induces anisotropic 
behavior in an originally isotropic alloy. 

(2) Because cast components are subjected to a cooling sequence designed 
to produce an anisotropic (single crystal or transversely isotropic 
alloy) . 

(3) Because a material can exhibit elastic anisotropy due to prior non- 
linear deformation. 

The application of the boundary element method to such materials requires a 
number of modifications. In particular, the following items require attention: 

(1) The point load solution for the (generally oriented) anisotropic ma- 
terial must be available analytically and calculable numerically. 

(2) Modifications must be made to the surface stress calculation. 

(3) Mew particular integrals for centrifugal body forces are required. 

Substantial progress was made during Task IIC in extending BEST3D to include 
anisotropic materials. 
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It has long been known that the point load solution for a general anisotropic 
material can be represented In terms of a line Integral, 

U.. (x - i) = — J (4) ds (3.4-16) 

1J “ Sir 2 | x-y| y 


where , 


= C ijkm ^k ^m 


C--b m = elastic constants 
ijkm 


and the integration path lies on the unit sphere and is perpendicular to the 
vector (x-y). 


In an earlier program (Reference 8), a technique was developed to tabulate the 
point load functions for directionally solidified (transversely isotropic) and 
eutectic (single crystal) materials. The tabulated results were then numeri- 
cally interpolated within a boundary element method code whenever kernel func- 
tion values were required. It was found that this method allowed accurate 
solution of stress analysis problems for both material types. The numerical 
evaluation of the kernel function was, however, considerably more expensive 
than the evaluation of the isotropic Kelvin solution, and led to an increase 
in solution time of about 100 percent. 


It is not known for exactly what forms of anisotropy equation (3.4-16) has a 
closed form solution. A condition for the existence of such a solution based 
on the root structure of a sixth degree polynomial involving the elastic con- 
stants has been conjectured (Reference 9), but not yet published. At present 
the only known closed form solutions are those for the isotropic (Kelvin solu- 
tion) and transversely isotropic cases. 


The transversely isotropic solution (Reference 10) was derived using stress 
functions involving complex variables. The original formulation was reviewed 
and considerably simplified before incorporation in BEST3D. In particular, the 
general form of Reference 10 includes within it a number of special cases. 
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only one of which is relevant to directionally solidified alloys presently in 
use in gas turbine engines. This fact was used to reduce significantly the 
amount of calculation required. 

The case of the single crystal material must still be handled using an approx- 
imate representation of equation (3.4-16). New methods for this approximation, 
requiring less storage and calculation than the technique reported in Refer- 
ence 8 are being developed as part of Tasks I VC and VC. 

In Task IC the effect of a centrifugal load on an isotropic material was eval- 
uated using a particular solution of the equilibrium equations. The particular 
solution has now been generalized for an arbitrary anisotropic material. If a 
body rotates about an axis passing through an origin, then the equilibrium 
equations are: 


r. . . 

iJ.J 



(3.4-17) 


where 

2 2 

_f = ( pu x, pcj y , 0) 
p = density 

u = speed of rotation. 

For simplicity in the presentation, the axis of rotation has been taken to be 

the z-axis. The problem can always be transformed to this case within the 

code. Such a transformation requires a corresponding transformation of the 

C.. matrix, which can then become a full, symmetric matrix, even for trans- 
1 J 

versely isotropic or single crystal materials. 

By analogy with the known particular solution for the isotropic case, the dis- 
placements of the anisotropic particular solution are assumed to have cubic 
variation in the Cartesian coordinates, for example: 

U 1 = A l x l 3 + A 2 x 1 2x 2 2 + A 3 x 1 x 3 2 (3.4-18) 
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where the coefficients are to be determined. Differentiation of (3.4-18) 
allows the determination of the strains (quadratic functions of the coordi- 
nates). The stresses are then determined, using the transformed Hooke's law, 
and a further differentiation allows the evaluation of the left-hand side of 
(3.4-17) in terms of the undetermined coefficients. Equation (3.4-17) is lin- 
ear in the Cartesian coordinates, and must be satisfied for an arbitrary 
choice of coordinates. This condition leads to a set of nine linear equations 
in the nine undetermined coefficients. This determination of coefficients is 
done only once for each subregion within a multiregion structure. Once the co- 
efficients in (3.4-17) have been determined, the anisotropic particular solu- 
tion is used in exactly the same manner as the isotropic solution. 

3. 4. 3. 3 Transient and Dynamic Stress Analysis 

3. 4. 3. 3.1 Review 

One of the goals of the advanced formulation development portion of the in- 
elastic methods program is the development of a three-dimensional boundary 
element capability for transient and dynamic analysis. Several different prob- 
lem types are of interest in the dynamic and transient analysis of gas turbine 
engine hot section components. The primary areas of interest are: 

(1) Determination of natural frequencies, and corresponding mode shapes, 
for geometrically complex structures. 

(2) Evaluation of the response of a structure to a periodic loading (par- 
ticularly near a natural frequency). 

(3) Determination of the time history of the response to a transient, 
nonperiodic loading. 

The first two problem types are related primarily to the avoidance of struc- 
tural problems due to forced vibration. The final type of analysis is normally 
used to predict the effect of unusual, and potentially very damaging, events 
such as impact by foreign objects. 
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During Task IC, the problem of dynamic/transient response was attacked by 
applying the boundary element method (BEM) to a transformed (Laplace or 
Fourier) form of the governing differential equations. Since a point load 
solution exists for the transformed equations, the overall framework of the 
static elastic solution could also be employed for the solution of dynamic/ 
transient problems. The version of BEST3D delivered at the end of Task IC in- 
cluded the capability for this transform domain solution. The proper execution 
and accuracy of this approach was demonstrated using a number of test prob- 
lems. The analytical and numerical development of this capability are dis- 
cussed in considerable detail in the First Annual Status Report (NASA 
CR-174700) . 

While the transform approach is capable of providing accurate solutions to 
dynamic/transient problems, it suffers from serious defects as a practical 
analysis method for complex structures. The most serious of these are: 

(1) Since the Laplace/Fourier transform casts the entire problem in the 
complex domain, the storage and computer time requirements are con- 
siderably increased. Typically, the solution for a single value of 
the transform parameter will cost two to four times a single static 
analysis for the same boundary mesh. 

(2) Any problem with nonperiodic loads must be solved by taking trans- 
forms of the loads, solving the transformed BEM problem for multiple 
values of the transform parameter, and then numerically inverting to 
reconstruct the time domain solution. Since the transform parameter 
is embedded nonli nearly in the point load solution, the BEM solution 
for each transform parameter value requires a complete reconstruction 
of the system equations, leading to extremely long computing times. 

(3) Since the frequency is involved nonlinearly, natural frequencies can 
be determined only through some form of a search procedure, a very 
expensive technique. 
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It was clear that, while the transform/compl ex variable method did solve the 
mathematical problem, new techniques would be required for a practical exten- 
sion of the BEM method to dynamic and transient problems. Very significant 
progress has been made In this area during Task IIC. In particular: 

(1) A real variable technique for the calculation of natural frequencies 
and mode shapes, originally developed for single region, two-dimen- 
sional analysis has been significantly improved, and tested in a two- 
dimensional BEM program. An extension of the method to multi region 
three-dimensional problems has been developed. 

(2) A real variable, time embedded formulation for the transient elastic 
problem has been developed and implemented in BEST3D. 

These developments are discussed in the following sections. 


3. 4. 3. 3. 2 Calculations of Natural Frequencies and Mode Shapes 


The governing differential equation is 


a 2 u. 

^ + ^ TxTajcT + 

l J 


3 2 U. 
J 


2 

+ pu) U . 

3X-3X. M J 


= 0 


which can be written in operator notation as 


l(u.) + pu u. = 0 
J J 


(3.4-19) 


The solution can be formally represented as the sum of a complementary func- 
tion w. satisfying 

J 


L(Wj) = 0 


(3.4-20) 


and a particular integral Vj satisfying 


L(v.) + pu u. = 0 
J J 


(3.4-21) 
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This equation, however, still contains the unknown displacement field 
which must be eliminated by a suitable approximation. Representing u^. as 



Uj(x) 


s c 4k (x,{") fUJ 1 ") 
m=l JK K 


(3.4-22) 


where Cj k is a function (like a kernel function) and d k are a family of 
known functions, equation (3.4-21) can be modified to give 

L(v.) + pu 2 2 c..(xj m ) d,(! m ) = 0 (3.4-23) 

3 m=1 JK k 


The function Cj k (x,l m ) is then chosen as 


c jk (x,l m ) = A « jk (1 - r m /A) (3.4-24) 

where A is a constant and r m is the distance between a boundary point and a 
set of reference points m on the boundary or in the interior. 

A particular integral can then be obtained in the form 

v . ( x ) = pu 2 2 D.JxJ m ) d(£ m ) (3.4-25) 

J m=l JK 

to satisfy equation (3.4-23). For the boundary points, the particular integral 
can be expressed as 


v = pu 2 D d (3.4-26) 

Once the displacement field due to the particular integral is known, the sur- 
face traction due to this displacement field can be expressed as 

i.V 2 t A 
t = pu 10 
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The final system of equations for the eigenvalue problem can be then con- 
structed using the particular solution (3.4-26) and the standard boundary ele- 
ment formulation as: 

Gt - Fu = pu 2 [G T * - F D d] = pu 2 [G T - F D] 6 (3.4-27) 

£ can be eliminated from equation (3.4-27) by invoking the matrix form of 
equation (3.4-22), i.e., 

u = C d or d = C" 1 u (3.4-28) 

Substituting (3.4-28) in (3.4-27), leads to: 

Gt - Fu = pu> 2 [GT - FD] C* 1 u = pu 2 M u (3.4-29) 

By incorporating the known and unknown boundary values, equation (3.4-29) can 
be recast into: 

Ax + po 2 Bx = 0 (3.4-30) 

where x is the vector of unknown boundary displacements and tractions, and A 
and B are two nonsymmetric fully populated matrices. It should be noted that 
the matrix A is identical to the system matrix of an ordinary elastic problem. 

3. 4. 3. 3. 3 Time Embedded Transient Dynamic Formulation 

The boundary integral formulation at boundary point £ for a general transient 
elastodynamic problem is given by 

[«i j - Cjj] Uj({,t) * [G^j*tj - F^j*u,j] ds(x) (3.4-31) 


where 


t 


G ij n i = l 6 i0 


(x,t;4,T) ^(x.t) dr 



★ 


u i 
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are Riemann convolution integrals where the body has been assumed to have been 
excited from rest. The functions and F^, which represent the dis- 
placements and surface tractions at a point (x,t) due to a unit force vector 
acting at a point | at a time t within a three-dimensional solid of infinite 
extent, are analytically too complex to reproduce here (see Banerjee and 
Butterfield) . 

Equation (3.4-31) represents an exact formulation involving integrations over 
the surface as well as the time history. Equation (3.4-31) is an implicit time 
domain formulation since the displacements at a time t are calculated taking 
account of the history of surface tractions and displacements up to and in- 
cluding the time t. Therefore, if grossly simplified assumptions are not made 
in the time variations of quantities, the stability problems encountered in 
the finite difference and the finite element methods should not arise. 


If the response at a time t is to be determined and the time domain has been 
represented by N nodes (with the node 1 at time t = 0 and the node N at time 
t), equation (3.4-31) can be rewritten as 
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where [***] denotes [G^ t^ - F^ u^]. 
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It is of interest to note that equation (3.4-32), like equation (3.4-31), is 
still an exact formulation of the problem since no approximation has yet been 
Introduced. 

By integrating analytically over the time Intervals and introducing a spatial 
discretization in the usual manner, (3.4-32) can be expressed as 

Ax( t) = By( t) - 2 [A N x N - B N y N ] (3.4-33) 

where the terms involved in the indicated sum include the effects of past 
dynamic excitation history on the boundary (i.e., from times t = 0 to t = 

Vl** 

If the time steps are kept constant, the left-hand side of the system equation 
needs to be generated only once, while the right-hand side is newly calculated 
for each time step, up to a maximum time. After this maximum time is reached 
no further integration Is required, no matter how many time steps are taken. 
The maximum time can be calculated based on the wave speeds of the material 
and the dimensions of the structure. Equation (3.4-33) can be written formally 
for each generic modeling region (GMR) as: 

Ax( t) = b(t) (3.4-34) 

which can be assembled for each GMR and solved for a multi -GMR problem in the 
usual manner. 

3.4.4 Computer Program Development 
3.4.4. 1 Introduction 

The computer program BEST3D was developed during Task IC to provide a tool for 
applying the three-dimensional boundary element method (BEM) to the elastic, 
inelastic and dynamic structural analysis of gas turbine engine hot section 
components. 
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The program* described in detail in the First Annual Status Report (NASA CR- 
174700), was designed to accommodate structures with very general geometry and 
loading. Further, it was clear that additional capabilities would be developed 
for BEST3D over a period of several years. For this reason the program was de- 
signed so that the anticipated capabilities could be incorporated within the 
original framework and data structure of BEST3D, without requiring major re- 
coding of already existing capabilities. 

During Task IIC, a number of new and enhanced capabilities were developed and 
installed in BEST3D. The basic structure of the program remains, however, very 
similar to that of the original version. The major changes and additions made 
to the code during Task IIC are described in the following sections. A full 
description of BEST3D may be found in the First Annual Status Report. 

Finally, it should be noted that development of BEST3D is carried out at Pratt 
a Whitney on an IBM 3084, in both batch and interactive modes, and on an 
HP9000 minicomputer at SUNY-B. The program is installed on the CRAY computer 
at NASA-Lewis Research Center. It can therefore be anticipated that BEST3D can 
be installed, with relatively little difficulty, on most commonly available 
computer systems. 

3. 4. 4. 2 Global Program Structure 

The major changes in the overall structure of BEST3D are: 

(1) The incorporation of a branch for the time embedded transient solu- 
tion. 

(2) Deletion of the complex variable option for transient analysis. The 
capability for complex variable analysis of forced response at a 
given frequency has been retained. 

(3) A branch has been designed for natural frequency/mode shape calcula- 
tions, and will be incorporated as part of Task I VC. 
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The overall structure of BEST3D at the end of Task IIC Is shown In Figure 
3.4-2. The program consists of a common Input section, followed by three 
branches, for static, forced response and transient analysis. The static 
analysis branch Is the model for the entire code, since the other branches 
largely employ generalized forms of the same algorithms used In the static 
analysis. The branch designed for natural frequency /mode shape calculation Is 
actually part of the static analysis loop. 

The transient and natural frequency /mode shape branches will be discussed 
separately. Other modifications, not visible at the level of the overall pro- 
gram structure (particularly in the plasticity analysis), will also be de- 
scribed. 

3. 4. 4. 3 Program Input 

The input to BEST3D is essentially unchanged. The discussion contained In the 
First Annual Status Report (NASA CR-174700) remains accurate. In particular, 
the library of surface elements has not been changed during Task IIC. The In- 
put changes made during Task IIC Include: 

(1) Very minor changes to provide the time step definition required for 
the time embedded transient analysis. 

(2) Simplification of the input used to describe nonlinear material 
models. 

(3) Inclusion of flags and (when required) additional input for embedded 
discontinuities, hot spots and rectilinear volume cells. 

Program input is fully described in the BEST3D Users' Manual. 
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figure 3.4-2 BEST3D - Overall Program Structure 





















3. 4. 4. 4 Surface Integral Calculation 


The most time consuming portion of an elastic or elastodynamlc analysis Is the 
numerical evaluation of surface Integrals. BEST3D uses product Gaussian quad- 
rature rules to calculate these Integrals. Considerable Improvement has been 
made during Task IIC in the efficiency and accuracy of these calculations. In 
addition, a fundamental difference between surface integration in the static 
and in the transient analyses has been clarified. 

3. 4. 4. 4.1 Integration in the Static Case 

In the original version of BEST3D, 4x4 product Gauss rules were used for all 
of the nonsingular surface integral calculations. Error estimates based on the 
one-dimensional estimates of Stroud and Secrest (Reference 11) were used to 
estimate truncation error in the integration. In the event that the requested 
error tolerance could not be met over the full element, a weighted rectangular 
subdivision of the element was used. During Task IIC, two things became clear. 

(1) A very large proportion of the time spent in numerical surface Inte- 
gration was used in 'near singular' Integrations, that Is, cases in 
which the source point was very close to, but not on, the element 
being integrated. 

(2) The (analytically derived) error estimates could sometimes be anti- 
conservative for low integration orders. They also tended to under- 
estimate the benefit of increased integration order. 


Based on these observations, studies were begun to improve the error estimates 
for the Gaussian quadrature. The objectives were first to develop more uni- 
formly accurate error estimates and then to employ these new estimates to re- 
duce the computational effort required to carry out the numerical surface In- 
tegration. 
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The new error estimates were developed by: 


(1) Noting that the error in the two-dimensional integration can be ex- 
pressed in terms of appropriate one-dimensional error estimates. 

(2) Carrying out the integration of traction kernel-shape function pro- 
ducts over a variety of randomly oriented source point-line segment 
pairs. 

(3) Fitting the results of these numerical experiments to a functional 
form based on the analytical error estimates and using statistical 
techniques to adjust the error equation so that the probability of 
exceeding the input error tolerance is less than 0.001. 

Comparison of the new and old error estimates led to several interesting con- 
clusions: 

(1) Higher order Gauss rules were found to be considerably more effective 
than previously thought. As a result, the library of integration 
rules in BEST3D was expanded to include orders two through thirteen, 
inclusive. This allows many more source point-element pairs to be in- 
tegrated without subdivision. 

(2) The ability to estimate error more accurately has allowed creation of 
a more effective algorithm for element subdivision. The original sub- 
division algorithm always performed a rectangular subdivision, using 
a 4 x 4 Gauss rule on each subelement. The new algorithm, in the most 
general case, uses a single central subelement surrounded by rings, 
each containing at most four subelements (Figure 3.4-3). The order of 
integration on each subelement, in each direction, is independently 
set to meet the requested tolerance. 
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Figure 3.4-3 Subdivision Strategy for Nonsingular Surface Integration 

The incorporation of the revised error estimate and subdivision algorithm in 
BEST3D has led to a reduction in computing time averaging 35 percent for elas- 
tic analyses, relative to the original version of the code. In addition, solu- 
tion accuracy is now much more closely related to the input error tolerance. 

3. 4. 4. 4. 2 Surface Integration in the Transient Algorithm 

The development of the embedded time algorithm for transient elastic analysis 
during Task IIC has led to a greatly improved understanding of fundamental nu- 
merical differences between the static and dynamic cases. In the static case 
it was found that improvement in integration accuracy and efficiency could be 
obtained by incorporating higher order Gaussian quadrature rules in BEST3D. It 
was originally anticipated that a similar improvement would result from incor- 
porating higher order rules in the transient branch of the code. It was found, 
however, that the introduction of the higher order rules led to greatly in- 
creased run times with little improvement in accuracy. 
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Further study clarified the fundamental difference involved. In the static 
case the integrands in all of the nonsingular surface integrals are infinitely 
differentiable and solution accuracy can, therefore, always be improved by the 
use of increased integration order. In the time embedded case, however, the 
point load solutions are only continuous. Physically this corresponds to the 
fact that the disturbance due to an impulse applied to a given time and spa- 
tial location is nonzero, at some later time, in only a finite portion of 
space. This means that the kernel function may be nonzero over only part of a 
given surface element. While the integrand is infinitely differentiable within 
both the zero and nonzero regions considered separately, its overall contin- 
uity over the entire element is only C°. The use of higher order quadrature 
rules is, therefore, of little use in improving accuracy. 

Based on these observations, a revised integration strategy was adopted for 
the transient branch of BEST3D. The surface elements are subdivided into a 
relatively large number of subelements, and relatively low (usually 2nd or 
3rd) order quadrature rules are used over each subelement. This has led to 
much improved accuracy in the transient analysis. 

3. 4. 4. 5 Volume Integrals 

No additional changes were introduced to the volume integration scheme de- 
scribed in First Annual Status Report (NASA CR-174700). 

3. 4. 4. 6 System Matrix Assembly for Multi region Plasticity 

For the single region inelastic problem the system equation can be expressed 
as: 


Ax = b 

where b contains the total effects of the boundary loading as well as the non- 
linearity. Since this final system matrix is algebraically identical to that 
for a single region elastic problem, the multi-GMR (generic modeling region) 
plasticity problem can use the same code for the assembly of the multi-GMR 
system equations. 
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3.4.4. 7 System Equation Solution 


The solution of the system equations Is essentially unchanged from the Initial 
version of BEST3D. The only major change was provision of a capability for 
more efficient solution process for problems Involving multiple loads. This 
was primarily required for efficiency In inelastic analyses. 

3. 4. 4. 8 Inelastic Solution Process for a Multi-GMR (Generic Modeling Region) 
Probl em 

The inelastic solution algorithm has been developed in such a manner that the 
existing coding for the multi region elastic problems is utilized without modi- 
fication. This essentially requires that the system matrix A and the right- 
hand side b be individually formed for each GMR and supplied to the solver for 
assembly and solution. The process of determining b of course requires the 
initial volume Integration of a number of algebraically complex kernels as 
well as the definition of the Inelastic states via the various constitutive 
models. This requires multiple solutions of the same system equations to 
satisfy the state dependent constitutive equations at all boundary and In- 
terior nodal points of cells. 

3.4.5 Val idation/Veri fi cation 

3.4.5. 1 Summary 

During Task IIC, significant new capabilities were added to BEST3D. The vali- 
dation and evaluation of these new capabilities and the continuing verifica- 
tion of BEST3D is discussed in the subsections which follow. Attention In this 
report is directed particularly to the validation and verification of the non- 
linear analysis capability and to the evaluation of the time embedded dynamic 
formulation. 
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A very large set of test cases was constructed during Task IC to validate 
BEST3D's ability to represent properly complex geometries and boundary condi- 
tions and to carry out accurately elastic stress analysis. Many of these test 
cases were executed again during Task I IC in order to verify the continuing 
correct operation of already existing capabilities. These test cases were dis- 
cussed in the First Annual Status Report (NASA CR- 174700), and are not dis- 
cussed further in this report unless they have been modified to demonstrate a 
newly implemented capability. 

3. 4. 5. 2 Validation of Elastic Capabilities 

Two major modifications were made to the elastic capabilities of BEST3D during 
Task IIC. First, the closed form point load solution for a transversely iso- 
tropic material was recast in a computationally efficient form and coded. 
Additionally, the required modifications to allow calculation of elastic 
stresses and strains at boundary nodes and to account for centrifugal ly loaded 
anisotropic materials were carried out. Second, a new error estimate for non- 
singular surface integration and a more efficient subdivision algorithm based 
on this estimate were coded and incorporated in BEST3D. The validation of 
these capabilities is discussed below. 

Transversely Isotropic Cube in Tension 

A unit cube was modeled using six linear surface elements, one for each face 
of the cube. The beam was loaded along the x-axis, parallel to one of the 
secondary material axes. The axis of isotropy was aligned with the z-axis. The 
displacements and stresses from the BEST3D analysis are compared to the exact 
solution in the table below. 


BEST3D Exact 


U x - elongation 
U - Poisson contraction 
U z - Poisson contraction 

a x 

All other stresses 


.0010003 in. 
-.00014321 
-.0050946 
21201 psi 
0 


.0010003 in. 
-.00014318 
-.00050942 
21200 psi 
<1 
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It Is clear that the boundary element method (BEM) results are essentially 
exact, validating the correct operation of the anisotropic capability. The 
analysis time for the transversely Isotropic case Is approximately 22 percent 
higher than that for the isotropic case. This is due to the fact that the 
transversely isotropic point load solution is computationally much more com- 
plex than the isotropic solution. 

The same geometry was also analyzed using two GMRs (generic modeling regions) 
in order to demonstrate the correct operation of the multi -GMR capability for 
anisotropic materials. 

Val i dati on of Surf ace Integration Algorithm 

The new surface integration algorithm was evaluated using the slab shown in 
Figure 3.4-4. Six boundary elements were used to model the slab. Although geo- 
metrically simple, the problem fully exercises the Integration algorithm be- 
cause of the large variation in the ratio of element size to the source point 
- element distance. 



Figure 3.4-4 Test Geometry for Surface Integration Algorithm 
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The variation of solution error (measured in terms of the reactions at the re- 
strained end of the slab) with the allowable relative error in the numerical 
integration is shown in figure 3.4-5. It is clear that the solution error is 
very directly related to the integration tolerance specified, which was not 
the case with the estimates used in the original version of BEST3D. As shown 
in Figure 3.4-6, it is also clear that the computational effort required to 
achieve a given level of solution accuracy grows very rapidly if very small 
tolerances are imposed. For this reason a default tolerance was set which 
attempts to optimize the relationship between accuracy and cost. 



Figure 3.4-5 Solution Error Vs. Integration Tolerance 
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COMPUTATIONAL EFFORT VS INTEGRATION TOLERANCE 



INTEGRATION TOLERANCE 

Figure 3.4-6 Computational Effort Vs. Integration Tolerance 


3. 4. 5. 3 Multi-region Nonlinear Analysis 

The multi -GMR nonlinear capability of BEST3D has been evaluated using an ex- 
tensive set of analyses of pressurized thick shells and disks. The results of 
these very demanding analyses, involving large-scale plasticity, were used 
both to validate the multi -GMR plasticity capability and to develop an under- 
standing of proper volume cell definition. The problems were run using an 
elastic-perfectly plastic material model and plane strain boundary conditions 
in order to allow comparison to an analytical solution. 
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Thick Cylinder 


The first problem discussed is that of a thick cylinder (outer radi us/inner 
radius = 2) subject to steadily increasing internal pressure. A 22.5° sector 
of the cylinder was modeled (Figure 3.4-7). Two boundary and volume discreti- 
zations, employing quadratic variation, were used. For both meshes the analy- 
sis was carried to collapse, that is, to the load at which BEST3D could no 
longer converge to a solution. 



MESH 1 

to DO UNOAftV CLEMENTS 
2 VOLUME CELLS 



MESH 2 

18 BOUNOAAV CLEMENTS 
A VOLUME CELLS 


Figure 3.4-7 BEST3D Models for Thick Cylinder 


The normalized load - deflection response of the cylinder is shown in Figure 

3.4- 8. Even Mesh 1, with only two volume cells gives a good representation of 
the solution, while Mesh 2 converges and gives an accurate displacement value 
up to 92 percent of the theoretical collapse load. The BEST3D results for the 
stress variation in the cylinder are also very accurate, as shown in Figure 

3.4- 9. 
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The BEST3D analyses for this problem were carried out using both one and two 
region models (with the interface located at a radius of 1.5a). The results 
were identical. It should also be noted that, in these analyses, the boundary 
elements and volume cells were equally spaced in the radial direction. Experi- 
ence in other problems indicates that convergence even closer to the theoreti- 
cal collapse load can be achieved by weighting the volume cells towards the 
inner radius of the cylinder. 


The second problem is that of a disk (outer radi us/inner radius = 5) subject 
to pressure at the inner radius. This geometry is typical of that encountered 
in high-pressure turbine disks. Experience with the cylinder analyses, dis- 
cussed above, indicated that at least four volume cells would be required for 
this analysis. It was also expected that an appropriate weighting of the cell 
discretization would prove desirable. The three meshes used (Figure 3.4-10) 
all consist of twenty quadratic boundary elements (including the interface 
elements) and four volume cells. 



UESH 3 


Figure 3.4-10 BEST 3D Models for Inelastic Analysis of Disk 
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As anticipated, the weighting of the volume cells toward the Inner radius 
leads to Increasingly good agreement of the BEST3D results with the analytical 
solution (Figure 3.4-11). It Is clear that more volume cells would be required 
to carry the analysis to collapse. In practice this would not be needed, since 
plastic deformation Is normally limited to the Inner 20 percent to 30 percent 
of a disk. 

It is clear (Figure 3.4-12) that all three meshes provide an accurate elastic 
solution. It is, therefore, of considerable Interest to study the optimization 
of the volume discretization for a fixed boundary model. This study is now 
being conducted. 

The analysis of the disk was repeated using a single region with cyclic bound- 
ary conditions. The results were unchanged from those discussed above. 



Figure 3.4-11 Load - Deflection Response for Disk 
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RADIUS/INNER RADIUS 

figure 3.4-12 Elastic Hoop Stress Vs. Radius 


3. 4. 5. 4 Application of Exact Initial Strain Solution 

The exact solution for the stress and displacement due to a uniform state of 
initial strain (or stress) in a rectangular parallel piped (discussed in Sec- 
tion 3) has been exploited in three different areas. A test case for one of 
these applications is discussed below. 




The capability to include the effects of embedded holes without explicit mo- 
deling has been incorporated in BEST3D. Up to ten such discontinuities may be 
included in each generic modeling region (GMR) . In order to evaluate the ef- 
fectiveness of this technique, for both elastic and plastic analysis, a test 
problem was defined in which a slot with rounded ends Is embedded in a tension 
strip (Figure 3.4-13). In the BEST3D analysis the slot is represented as a 
rectilinear cell, without explicit modeling. 

The BEST3D elastic results, for the stress variation between the free surface 
and the slot surface, are compared with a plane strain solution (using modeled 
slots) in Figure 3.4-14. The stresses show good agreement except near the slot 
surface, where the actual surface geometry of the slot becomes important. 



Figure 3.4-13 Constant Thickness Plate With Through Thickness Slot 
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Figure 3.4-14 Stress Variation in Plate With Slot 


3. 4. 5. 5 Embedded Time Dynamic Analysis 
Explosion in a Spherical Cavity 

The time embedded transient solution was exhaustively tested using the problem 
of a spherical cavity in an elastic space, subjected to a suddenly applied 
pressure load. Three boundary meshes (Figure 3.4-15) were used for the solu- 
tion of this problem. 

This problem was chosen for detailed study because the existence of an exact 
solution allows precise evaluation of the accuracy of the numerical solution 
and, in particular, the effects of time step size and boundary mesh refine- 
ment. Figure 3.4-16 compares the results of a BEST3D analysis, using Mesh 1 
and a time step of 0.00035 second, with the analytical solution. The boundary 
element method (BEM) results show good qualitative agreement with the exact 
solution, but underestimate the peak deflection by about 12 percent, due pri- 
marily to an overly crude boundary mesh. This is demonstrated in Figure 3.4- 
17, where the dynamic magnification factors (transient/steady state displace- 
ment) for the exact and numerical results are compared and quite good agree- 
ment is apparent. 
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EXPLOSION IN SPHERICAL CAUITY - P1ESH 1 



Figure 3.4-16 Explosion in Spherical Cavity - Mesh 1 



Figure 3.4-17 Explosion in Spherical Cavity - Dynamic Magnification 




The question of sensitivity to time step size Is always of Importance In 
transient analysis. Mesh 1 was also used to examine this question by calcu- 
lating the solution for a variety of step sizes (0.0002 second to 0.0007 
second). The results show almost no sensitivity to step size. 

Based on the experience with Mesh 1, it was concluded that a more refined sur- 
face mesh, capable of giving a more accurate result for the static problem, 
would also be suitable for the transient analysis. Mesh 2 and Mesh 3 were used 
to confirm this, although only the results for Mesh 3 are presented in the re- 
port. The dynamic magnification factors obtained with Mesh 3 are slightly im- 
proved relative to Mesh 1. The absolute displacement - time response, however, 
is dramatically improved (Figure 3.4-18). 



Figure 3.4-18 Explosion in Spherical Cavity - Mesh 3 
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3.4. 5.6 Calculation of Natural Frequencies and Mode Shapes 


A number of two-dimensional test problems were solved in order to validate the 
algorithm designed for natural frequency and mode shape calculations in 
BEST3D. One of these problems, which was also solved using the MHOST finite 
element code, is described below. 

The first four bending modes for a cantilever beam were calculated using both 
a two-dimensional BEM code and the MHOST finite element code. The beam has a 
length of 6.5 units and a square (lxl) cross section. The meshes used for the 
two analyses were very similar, with fourteen nodes along the beam length for 
the finite element model and thirteen for the BEM model. The calculated eigen- 
values agree well. 


Natural Frequencies (Hz) of Cantilever 


Mode 

2-D BEM 

MiOST 

1st bending 

0.368 

0.378 

2nd bending 

2.214 

2.188 

3rd bending 

5.591 

5.583 

4th bending 

9.986 

9.908 


Further, the mode shapes calculated using the two techniques are indistin- 
guishable. The first and fourth modes are shown in Figure 3.4-19. It should be 
noted that the fourth mode displays a nonzero slope at the fixed end. This is 
a real feature of the two-dimensional solution which is not present in a beam 
theory analysis, since it is suppressed by the fixed end boundary conditions 
normally used in beam theory. 
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Figure 3.4-19 Comparison of Finite Element and Boundary Element Mode Shapes 








3. 4. 5. 7 Inelastic Analysis of the Benchmark Notch Specimen 


The Benchmark Notch Specimen is a double edge notch specimen developed by 
General Electric/Louisiana State University (GE/LSU) under NASA-Lewis Contract 
NAS3-22522. A significant volume of well documented data was provided as part 
of the contract report (Reference 12). These data were used during Task IC to 
verify the elastic capabilities of BEST3D. During Task II, the inelastic data 
from this program have been used to verify the nonlinear capabilities of both 
BEST3D and MHOST. 

The specimen geometry is defined in Figure 3.4-20. Stress analysis was carried 
out for the gage section only, a procedure already known to be satisfactory. 
The specimen models used are shown in Figure 3.4-21 (for BEST3D) and Figure 
3.4-22 (for MHOST). In both cases it was found that selective mesh refinement 
and weighting of elements toward the notch were required in order to obtain 
acceptable accuracy with reasonable computing times. It should be especially 
noted that the MHOST models use eight node elements. This means that the MHOST 
“fine mesh" and "intermediate mesh 2" provide the same through- thickness dis- 
placement variation as the quadratic BEST3D elements. This refinement was 
needed to obtain accurate finite element results. 

In Figures 3.4-23 and 3.4-24, the calculated strains at the notch root are 
compared to a regression fit to the GE/LSU data. It is clear that, with an 
appropriate choice of mesh, both BEST3D and MHOST accurately predict the 
measured strains. 

In addition to the monotonic loading results discussed above, the cyclic 
plasticity model in BEST3D was used to predict specimen response over a more 
complex loading sequence (no load --> max load --> min load --> max load). The 
results of the cyclic calculation for Test 7 are shown in Figure 3.4-25. 
Again, excellent agreement with the GE/LSU data was obtained. 
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Figure 3.4-20 Double Edge Notch Specimen Used In Contract NAS3-22522 
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Figure 3.4-21 


Boundary Element Models for Benchmark Notch Analysis 



Figure 3.4-22 


Finite Element Models for 


Benchmark Notch Analysis 
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figure 3.4-23 Benchmark Notch - BEST3D Vs. Measured Strains 



figure 3.4-24 Benchmark Notch - MHOST Vs. Measured Strains 
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figure 3.4-25 Benchmark Notch - First Cycle Response 
3. 4. 5. 8 Component Analysis 

Further elastic analysis has been done for the turbine blade (Figure 3.4-26) 
originally discussed in the First Annual Status Report (NASA CR-174700). The 
blade was analyzed under centrifugal load using both BEST3D and MARC (three- 
dimensional finite element). The evaluation of the results was greatly sim- 
plified through the use of an in-house program providing (color) deflected 
shape and iso- stress plots. It is a relatively straightforward task to link 
BEST3D to such a tool . 

In Task IC the intent of the turbine blade analysis was primarily to verify 
the correct operation of the code for this highly complex problem. The intent 
in Task I IC was to verify BEST3D's ability to calculate correctly the peak 
stress in the blade and determine the degree, if any, of mesh refinement 
needed for this calculation. 
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figure 3.4-26 BEST3D Model of Cooled High-Pressure Turbine Blade 


A fully linear BEST3D analysis correctly predicted the peak stress location 
and gave a peak (maximum principal) stress of 146 ksl, significantly lower 
than the MARC result of 169 ksi. Experience with other problems Indicated that 
the use of limited quadratic variation near the peak stress location should 
lead to considerable improvement with a very limited increase in computing 
cost. 
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In order to test this supposition, the 8EST3D analysis was repeated using 
eighteen quadratic elements. The extent of quadratic variation is indicated in 
Figure 3.4-27. This led to an increase in the peak stress to 174 ksi, in ex- 
cellent agreement with the MARC result (which is calculated at an element in- 
tegration point somewhat below the surface). The increase in computer time was 
only 7 percent, i.e., from 890 to 952 central processing unit (CPU) seconds. 
By comparison the MARC analysis required over 2700 CPU seconds to obtain 
equivalent results. 



Figure 3.4-27 Extent of Quadratic Variation Used in BEST3D Turbine Blade 
Analysis 
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3.4.6 List of Symbols 



List of Symbols 
Referenced Within Section 3.4 


Symbol 

Description 

Page 

x j 

Cartesian coordinates 

132 

. o 

Initial stress rate 

133 

X,u 

Lame constants 

133 

u i 

Displacement rate 

133 

6 ij 

Kronecker delta 

133 

£ 

Young's modulus 

146 

a 

Coefficient of thermal expansion 

133 

T 

Temperature 

131 

f 

Body forces loading 

131 

S 

Boundary of body to be analyzed 

131 

S 

Point on boundary S 

132 

X 

Integration point 

132 

yi 

= x i " 1 

132 

r 

= y 

132 

V 

Interior of body to be analyzed 

131 

G U 

Kelvin solution 

131 

Tijk 

Stresses derived from G-jj 

131 

FlJ 

Tractions derived from T^j^ and 
surface normal 

131 

C 1j 

1/2 If S Is smooth at y, 

otherwise depends on surface geometry at 

131 
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List of Symbols 

Referenced Within Section 3.4 (continued) 


Symbol 

Description 

Page 

*i 

Tractions 

131 

°jk 

Stress rate 

131 

Dijk 

Sijk 

E ipjk 

Higher order kernels derived 
from Gjj by differentiation 
and use of Hooke's Law 

131 

• 

X 

Vector of all unknown freedoms 
(displacements and tractions) 

134 

• 

y 

Vector of all known freedoms 
(displacements and tractions) 

134 

p 

Mass density 

143 

u 

Frequency 

143 
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